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I. INTRODUCTION 


A. MOTIVATION 


The continual development of computer technology has enabled the expansion of 
intelligent control into the field of remotely operated underwater vehicles. With this 
advent, tetherless Autonomous Underwater Vehicles (AUVs), have arrived and are 
expected to be self-sufficient with respect to power and control. Driven by the need to 
lower the cost of operations, applications associated with these vehicles are both 
appealing and numerous. The potential uses for these vehicles include but are not limited 
to: scientific (oceanography, geology, geophysics,...[Curtin 1993, Smith 1994 and 
Pereira 1996]), environmental (waste disposal monitoring, wetland 
surveillance,...[Chase 1998]), commercial (oil and gas, submerged _ cables, 
harbors,...[Hartley 1991, Butler 1998, Kato 1998]), or military (mine warfare, tactical 
information gathering, smart weapons,...[Honegger 1996]). During a mission, an AUV 
is expected to carry sensors (side scan sonar, bathymetry, bottom profiler,), track to a 
certain planned trajectory (traveling from point A to point B, performing the systematic 
exploration of a zone,...), and even make on-line decisions allowing for mission 
reconfiguration, [Yavani 1996, Byrnes 1996]. 

The payload, which the vehicle carries, can place constraints on the motion of the 
vehicle and mission planning. With new emphasis on naval mine countermeasures and 
reconnaissance in very shallow waters (VSW) [ONR 1998], a typical mission may 
require a vehicle to follow an unknown and profile-varying bottom at a desired altitude, 
while maintaining a nominal constant forward speed, and avoiding obstacles when 
necessary, all while operating in a highly energetic environment. Since small AUVs are 
sensitive to wave induced motions [An 1996], AUVs in VSW will perform better if they 
derive knowledge, using their onboard sensors, about their environment and operating 
area. 

The question then is: How can environmental knowledge be obtained and used in 
a real-time control architecture in order to correctly and safely achieve an assigned 


mission? The integration of sensory data in closed-loop control of mobile robots or 
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vehicles has been a topic of recent research in both control architecture [Saridis 1983, 
Healey 1996, Marco 1996] and execution [Marco 1998]. Based on this work, two broad 
categories can be distinguished: 1) the behavioral approach and 2) the sensor based 
control approach. The first is too restrictive for general use, since it needs to specify the 
problem as a set of simple input/output relations for all possible situations. Among the 
second class, there are several methods proposed, including intelligent servo control 
[Marco 1996c], path following [Fryxell 1995] and collision avoidance [Williams 1990], 
and the task function approach [Santos 1995]. 

The task function approach allows a framework by which various control laws 
may be develop for specific mission tasks. The embedded control laws are then used 
based on sensor input decisions with a mission manager controlling the commands and 
task assignments. The goals of this research include the need to show that task 
assignments can be performed with less environmentally induced motions if the 
characteristics of the operating environment are known. This is particularly important 
because of the need to operate in environmentally energetic areas and because of the 
unpredictable nature of the VSW region. 

Classically, this problem falls into the broad category of servo control disturbance 
rejection. This problem is well known and has a rich history of study. In principle, a 
high gain feedback servo is known to track commands and reject disturbances (unwanted 
inputs). Requirements for loop stability and sensor noise rejection are competing and the 
subject of optimal control formulations. It is conjectured that knowledge or direct 
measurement/estimation of these disturbances improves control performance. 
Measurement/estimation of ocean disturbances is both difficult, as well as subject to the 
unpredictable nature of the ocean. For all the ongoing work in the offshore oil industry, 


wave disturbances are still not measured and used directly for control purposes. 


B. BACKGROUND 


Autonomous systems for small unmanned untethered underwater vehicles have 
several features that separate them from traditional marine control and guidance systems. 


The single most important feature is that it is desirable to control a small underwater 
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vehicle with relatively high velocities along or about two or more axes at the same time 


reliably. This leads to stronger coupling, larger nonlinearities, more state equations and 


more unknown parameters in the vehicle's equations of motion than what would be 


present with surface vessels. This single most important fact makes the control of small 


underwater vehicles relatively difficult. 


Other factors include but are not limited to: 


A small underwater vehicle may be controllable in all six degrees of freedom, 


A small underwater vehicle has a natural frequency in the same range as the 
environmental disturbances; 


Actuator dynamics are much faster on small underwater vehicles; 


Power for control and operations is limited by what the vehicle can carry 
onboard; 


Man-in-the-loop operations and human intervention for fault response are not 
possible, if controller problems develop; 


Small underwater vehicles, having a higher bandwidth, may compensate for 
the first order wave disturbances by means of feedforward and feedback 
control laws. 


There are numerous factors that must be considered in the design of a control 


system that is capable of being implemented in an operating vehicle. Some of the most 


important items are: 


The bandwidth is limited by the control system's sampling rate, computational 
time delays and sensor communication delays. The time delay factor is 
extremely important when the control law is dependent on information from a 
hydro-acoustic source. 


The control system's sensitivity to noise has to be considered, because it 
determines how good a sensor must be and how robust the control algorithm 
is with respect to noisy measurements. 


The required accuracy of the system is usually a function of the commanded 
tasks; hence it is not meaningful to discuss this factor before a task 
specification is made. 


e Stability is of the utmost importance. Depending on the type of system (SISO 
or MIMO) and the type of controller to be employed (linear or nonlinear), 
various techniques are available to rate these criteria. 


e Requirements to handle system failures, i.e., sensor failure, actuator failure 
and even computer failure, must be addressed. 


e Parameter variation is also an important aspect in robustness analysis of a 
control system. Parameter variations can be handled by a well-designed 
adaptive algorithm or by employing a control scheme that is robust to these 
variations. 


The first modern DP system was first introduced by [Balchen et al., 1976], where 
a Kalman filter approach for solving the wave noise-filtering problem was employed. 
Before the employment of this system, DP systems were mainly based on PID-controllers 
and matching notch filters which filtered out the wave "noise" on the sensors. Kalman 
filter based DP systems now dominate the market and the literature [Balchen 1976, 
Jenssen 1980, Saelid 1983, Fung 1983, Triantafyllou 1983, Fossen 1995, and Sorensen 
1996]. | 

The Kalman filter based DP systems consist of one high and one low-frequency 
model where the control action is based on the low frequency model. The high frequency 
model is used to filter out the relatively high-frequency, first-order, wave noise that 
appears on the sensors. 

This filtering is not necessary for small underwater vehicles, since in theory, 
compensation for the first order wave disturbance’s is possible. For small vehicles, the 
DP case becomes a special case of the general control problem where desired trajectory 
velocities are zero. Although there exist many available DP systems for surface vessels, 
the Simrad-Albatross as one example, and even several ROV DP systems, such as the 
[Marquest 1991] system and the systems developed and deployed on the Norwegian 
Experimental Remotely Operated Vehicle (NEROV) by [Fossen and Satagun 1991], 
there currently ae NO KNOWN AUTOMATIC DP systems for AUVs. 


Cc; OBJECTIVES 


The main objective of this research is to show through modeling, simulation and 
experimental validation that intervention tasks performed by intelligent underwater 
robots are improved by their ability to gather, learn and use information about their 
working environment. 


Specifically this dissertation addresses the following topics: 
e Generation and verification of mathematical models; 


e Measurement and use of wave disturbance information for control 
compensation; and 


e Implementation and verification of real-time embedded control processes. 


D. CONTRIBUTIONS 


This dissertation contains both theoretical and experimental contributions to the 
field of applied underwater vehicle control. The theoretical contributions are partly 
found in the modeling chapter (Chapter I), the disturbance rejection chapter (Chapter V), 
and the wave direction estimation chapter (Chapter VII). The experimental contributions 
(Chapter VIII) are from work carried out during operational testing of the NPS Phoenix 
AUV. The Phoenix and its subsystems, including system upgrades necessary to perform 
this research, are described in Appendix C and in [Marco 1996]. 

The contributions in this dissertation are described below: 

e Chapter II provides a new generalized approach to the modeling of small 
underwater vehicles subject to shallow water wave and current effects. Using 
appropriate modeling formulations, as opposed to adding white or colored 
noise, random disturbances, or "RAO" based motions, the fluid effects are 
incorporated directly into the equations of motion. In addition, this 
formulation provides the ability to construct a generalized distributed 


simulation capability for the evaluation of underwater vehicle systems in 
shallow water. (Generally useful to U.S. Navy tactical decision making). 


e It is proven through simulation (Chapter V), and verified by experimental 
validation (Chapter VIII), that measuring the water column velocities directly, 
wave induced disturbances can be substantially compensated by the vehicle 


control system. This technique eliminates the need to develop and incorporate 
sophisticated predictive disturbance models in the control system design. 


In Chapter VI, it is shown how small underwater vehicles, using direct fluid 
measurements, can obtain short-term wave magnitude, directionality and 
current estimates, thereby providing useful information in the area of tactical 
oceanography. It is also shown how a general seaway direction may be 
obtained for use in mission planning and control. 


E. DISSERTATION ORGANIZATION 


This dissertation is aimed at contributions in the area of control of intelligent 


underwater robots. Figure 1.1 represents the general control architecture for the Phoenix 


AUV. Using this figure as a roadmap, the chapters and appendices consider the 


following: 


Chapter IJ. "Mathematical Modeling Of Underwater Vehicles" derives the 
complete set of nonlinear equations of motion for a small underwater vehicle 
subject to shallow water waves. Kinematics, Newton's laws of angular and 
linear momentum, general hydrodynamics and external force modeling are 
discussed in detail. 


Chapter III. "Disturbance Analysis" describes how deterministic and 
stochastic disturbances can be modeled for use in the vehicle control 
architecture. Statistical description, state space representation and 
autoregressive (AR) modeling are used to illustrate these ideas. 


ChapterIV. "Parameter Identification" describes the theoretical foundation 
and experimental results associated with determining the parameters and 
coefficients used to model the NPS Phoenix AUV in the surge direction. 
Comparisons between vehicle experiments and simulations show the level of 
certainty associated with the identified parameters. 


Chapter V. "Disturbance Rejection Theory" provides a survey of the 
classical and modern approaches to disturbance rejection and compensation. 
This chapter describes how a nonlinear surge controller can be applied to an 
underwater vehicle design. It details three specific controller designs: a LOR 
design, a LOR with embedded disturbance model design, and a nonlinear 
sliding mode controller with disturbance compensation. 
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Figure 1.1 Phoenix AUV Control Architecture 


e Chapter VI. "Disturbance Compensation Controller (DCC)" is a discussion 
of the design and implementation of a combined nonlinear model-based 
extended Kalman filter with a SMC for dynamic positioning. The state 
estimator is critical to the controller performance since the surge controller 
requires states that are not measured with existing systems and sensors. This 
chapter compares actual signals with estimated signals, and shows that by 
properly designing the filter gains, the measurement noise levels are reduced 
and the controller bandwidth is increased. Also, critical to the point of the 
dissertation, ocean fluid particle velocities, the disturbance, are estimated. 


e Chapter VII. "Wave Directional Estimation From A Moving Platform" 
outlines the theory used to estimate wave directions from an AUV. It also 
describes the method by which a heading command, based on the dominant 
wave energy direction, is calculated. This heading command is used in the 
heading controller to properly orient the vehicle in the direction of the 
maximum energy, thereby reducing the drag force on the vehicle. 


e Chapter VIII. "Experimental Results and Validation" is an analysis of the 
experimental missions performed during this research. It describes the 
real-time implementation of the various processes (filters, controllers and data 
acquisition code) operating in the Phoenix. The chapter shows that the 
Disturbance Compensation Controller (DCC) is capable of dynamic station 


keeping in the presence of waves. It also displays and analyzes directional 
wave spectrum estimates obtained during AUVFEST '98. 


Chapter IX. "Conclusions and Recommendations" contains comments, 
conclusions and recommendations for future work. 


Appendix A. "Equations of Motion and Parameters for the NPS Phoenix 
AUV" provides the physical characteristics, the hydrodynamic coefficients 
and the equations of motion for the NPS Phoenix AUV. 


Appendix B. "Doppler Sensors" outlines the specifications and principles of 
operation of the SonTek® ADV and the RDI® WORKHORSE DVL. 


Appendix C. "The Phoenix AUV" describes in detail the vehicle more 
closely, including upgrades and acquisitions that allowed the Phoenix to 
transition from a test tank vehicle to open ocean operations. 


Appendix D. "AUVFEST '98" describes the NAVO sponsored underwater 
vehicle demonstration exercise that took place during November 1998 off the 
coast of Gulfport, MS. 


Il. MATHEMATICAL MODELING OF UNDERWATER VEHICLES 


A. INTRODUCTION 


The mathematical modeling of underwater vehicles can be found throughout the 
literature. The standard equations of motion for submarine simulation from the David 
Taylor Model Basin Naval Ship Research and Development Center (DTMSRDC) are 
described in detail by [Gertler 1967] and [Feldman 1979]. In [Kalske 1989] a survey of 
the motion dynamics of underwater vehicles, including ROVs and submarines, is given. 
A description of the linear and nonlinear motion dynamics of the University of New 
Hampshire Experimental Underwater Vehicle (UNH-EAVE) is found in 
[Humphreys 1982]. 

Many other papers and books have been written on this subject (see Yuh 1990, 
Healey 1992 and 1993, and Fossen 1994 for further examples). Each model has 
developed in complexity and accuracy, but each model has vehicle specific 
simplifications making it necessary to revisit this topic. Presently, there is no single 
model that combines all aspects, including environmental disturbances and their effects, 
into one generalized model for shallow water operations of small underwater vehicles. 

Accurate modeling allows for the development of precision control. Precise 
control is needed for many maneuvers and tasks, such as autonomous docking and 
recovery, target detection, identification and localization, as well as station-keeping. 

In this chapter the generalized six-degree of freedom (6DOF) equations of motion 
(EOM) for a small underwater vehicle operating in shallow water will be developed. As 
with all previous modeling work in this area, the underlying assumptions are that: 

1) The vehicle behaves as a rigid body; 

2) The earth's rotation is negligible as far as acceleration components of the 

center of mass are concerned; 

3) The primary forces that act on the vehicle have inertial, gravitational, 

hydrostatic and hydrodynamic origins, and 


4) The hydrodynamic coefficients or parameters are constant. 


The chapter will begin by outlining the coordinate frames and the kinematic and 
dynamic relationships used in modeling a vehicle operating in free space. Next a 
discussion of linear wave theory and basic hydrodynamics will be presented. This 
discussion will set the foundation for the various force and moment expressions 
representing the vehicle’s interaction with its fluid environment. The control forces, 
resulting from propellers and thrusters and from control surfaces or fins, that enable the 
vehicle to maneuver, will be then be detailed. With the hydrodynamic and control force 
and moment analysis complete, the full six degree of freedom equations of motion will be 
formed taking into account the necessary modifications for current and wave effects. 
While reduced order models representing flight control modes are detailed elsewhere, the 
chapter will conclude with a development and discussion of the one degree of freedom 
(1DOF) surge model. This model will be used in later chapters as the basis for a 


controller that will allow an AUV to station-keep in the presence of waves. 


B. COORDINATE SYSTEMS, POSITIONAL DEFINITIONS AND 
KINEMATICS 


For underwater vehicles that operate in three dimensional space and time, it is 
necessary to describe position/orientation and velocity/rotation rate of the vehicle by six 
independent coordinates or degrees of freedom. Typically these coordinates are chosen to 
correspond to the position and orientation and their time derivatives with respect to some 
set of mutually orthogonal coordinate axes fixed to an arbitrary origin which defines a 
reference frame. Likewise, the forces/moments acting on or produced by the vehicle can 
be referenced to a set of coordinate axes. In this dissertation, standard notation, 
[SNAME 1950], is used to describe the aforementioned 6 DOF quantities and is 
summarized in Table 2.1. 

Note that by convention for underwater vehicles, the positive x-direction is taken 
as forward, the positive y-direction is taken to the right, the positive z-direction 1s taken 
as down, and the right hand rule applies for angles. It is convenient to group the linear 
and angular body fixed velocities into a vector quantity x, where x = [u, v, w, P, q, ae 


and the global positions and Euler angles as a vector quantity z, where 
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z=([X, Y,Z, ¢, @ Wi’. The externally applied forces and moments are grouped into a 


vector quantity f, where f = [X;, Y, Zp Kp, Mp Nol’. 


Forces and Linear and Positions and 


Moments Angular Euler Angles 


Velocities 





Table 2.1 Standard Underwater Vehicle Notation 
1. Reference Frames 


As discussed earlier and outlined in Table 2.1, to properly describe or model the 
motion of a rigid body three independent positions and angles are required, and for 
convenience, three orthogonal coordinate frames are used. First, a global frame OXYZ, as 
shown in Figure 2.1, is defined with origin O, and a set of axes aligned with directions 
North, East and Down. This produces a right-hand reference frame with unit vectors I, J, 
and K. Ignoring the earth's rotation rate in comparison to the angular rates produced by 
the vehicle's motion, it can be said that the OXYZ coordinate frame is an inertial reference 
frame in which Newton's Laws of Motion will be valid. A vehicle's position in this 
global frame will have the vector components, ro: = [XI + YJ + ZK]. 

Secondly, a body fixed frame of reference Oxyz, with the origin O, and unit 
vectors 1, j, k, located on the vehicle centerline, moving and rotating with the vehicle is 
defined. The origin O’, will be the point about which all vehicle body force will be 
computed in later sections of this chapter. The vehicle's center of gravity (mass), CG, 


and center of buoyancy, CB, do not generally lie at the origin of the body fixed frame, nor 


1] 


are they collocated. The positional vectors of the CG and CB relative to the origin of the 
body fixed frame are Pg and Pp, respectively, and can be represented in component form 
as [xgi + ycj + zgk] and [xgi + yaj + zek]. The location of the center of mass (gravity) is 
important because Newton's Laws of Motion simplify when forces and moments on a 
body are applied to its center of mass. The center of buoyancy is determined by the 
shape of the submerged portion of the body, while the center of gravity is determined by 
the distribution of the weight. 

Lastly, a fluid frame is defined with origin O,. This frame is aligned parallel to 
the global reference frame but moves with the local fluid particles. Defining the fluid 
reference frame in this manner simplifies the hydrodynamic force modeling which will be 


discussed in later sections. 


GLOBAL FRAME 
O 
xX 
MW 
Z ro 
p 
O; G O ’ 
C DYDY-FIXED FRAME 
U; a p : 
y & 
Vr 
Wr 
FLUID FRAME 
(PARALLEL TO GLOBAL FRAME) Zz 


Figure 2.1 Coordinate Frames and Axes Convention 


De Euler Angles 


The transformation from one Cartesian coordinate system to another can be 
performed by means of three successive rotations in sequence. While there are several 
different methods to describe the attitude of a vehicle in the global reference frame, the 
most common method uses so-called Euler angles. This method uniquely defines the 
angular orientation of the vehicle reference frame relative to the global reference frame. 
There are several Euler angle conventions to describe these three successive rotations, see 
[Craig 1989] for a complete list. However, for the purpose of this dissertation, the so 
called "roll, pitch, yaw,’’ "X-Y-Z fixed angle" or "Tait-Bryant" convention will be used. 
One concern that may arise with the use of Euler angles is that a singularity exists when 
the elevation reaches 90 degrees. This limitation which can sometimes, although rarely, 
cause trouble in flight simulations and control computations may be overcome by the use 
of quaternions [Craig 1989], which introduce four rather than three variables to describe 
angular position. 
| For the "roll, pitch, yaw" convention, a forward transformation is performed by 
beginning with a vector quantity originally referenced in the global reference frame. 
Then, through a sequence of three rotations it is transformed into a frame that is assumed 
to be parallel] to and moving with the vehicle body coordinate frame. To start the 
transformation, begin by defining an azimuth rotation y, as a positive rotation about the 
global Z-axis. Next define a subsequent rotation 6, (positive up) about the new Y-axis, 
followed by a positive rotation ¢, about the new X-axis. The tmple rotational 
transformation in terms of these three angles, is then sufficient to describe the angular 
orientation of the vehicle at any time. 

It follows that any position vector, ro, in an original global reference frame given 
by ro = [Xo, Yo,Zo], will have different coordinates in a rotated frame when an azimuth 
rotation by angle y, is made about the global Z-axis. 

If the new position is defined by r; = [X;, Y;,Z;], it can be seen that there is a 
relation between the vector's coordinates in the new reference frame with those that it had 


in the old reference frame. It follows that, 
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X,=X,cosy+Y, siny (2a) 
Y,=-X,snw+t+Y, cosy, (22 
with Z;=Z,. This relation can be expressed in matrix form by the rotation matrix 
operation, 
ry =[Tyz Iro; (2.3) 
where the rotation matrix Ty,z, represents an orthogonal transformation. Multiplication 
of this rotation matrix with any vector, ro, will produce the components of the same vector 
in the rotated coordinate frame. Continuing with the series of rotations recut in a 
combined total rotational transformation, 
T1(9,6 w) = T() T(@) Ty). (2.4) 


In expanded notation Equation 2.4 takes the form 


cwc@ swc@ —s0 
T,(¢.0,y)= cwsso—swcgd sws0s6+cwcd cOso 
cws&o+swsd sysko-—cysd cco 


where c and s are abbreviations for cos and sin. It can be said that any positional vector 
in an original reference frame may be expressed in a rotated frame with coordinates given 
by the operation, 

rik = T1(9,8 W rux. (2.5) 


3. Kinematics 


Kinematics deal with the relationships of motion quantities regardless of the 
forces induced by their prescribed motions. Description of a body's position both 
translational and rotational, will need to be related to velocities, both translational and 
rotational, prior to extending the analysis to accelerations. The connection between 
translational velocity and the rate of change of translational position is straight forward. 


Define a global velocity as, 


P= ie (2.6) 


14 


This vector will have components that are different when seen in a body fixed frame. 
Selecting the linear components of the body fixed velocity vector, v = [u, v, w], these 
components, in global quantities are found, using the forward transformation defined in 


Equation 2.5, to be 


u se 
v1=7,(0,0,w) Y |. (2.7) 
w Z 


It is a simple reverse coordinate transformation from body fixed to global coordinates to 


see that, 
x Tu 
Y |=7, (9,6,y)| v |, (2.8) 
Z w 


where T,1(¢,0,y) is the transpose of Equation 2.4, since T,(¢,0,y) is orthogonal. This 
shows that the progression of a vehicle in the global frame clearly depends on its local 
velocity components and its attitude. | 

The connection between angular attitude and angular velocity is not quite so 
obvious. At first sight, it is tempting to define the instantaneous angular velocity of the 
vehicle simply as the rate of change of its angular position defined by the Euler angles. 
This is erroneous however, because the rotation @, was defined as a rotation about the 
intermediate frame after a rotation y, had been made. It should be noted that the order of 
multiplication in Equation 2.4 must be followed since the rotations do not commute, i.e., 
T (¢,0,w) #T(@,w,). Since the rotation do not commute, the set cannot form a vector 
space and therefore the derivatives cannot represent the body's angular velocity. 

Vehicle inertial angular rates are defined in terms of components that have 
angular velocities about the global axes. It is necessary to relate both Euler angles and 
their rates of change, to angular velocity components about the global axes and to their 
components lying along the body fixed axes in any attitude. The prime reason for this is 
that it is difficult to construct physical sensors to measure rates of change of Euler angles. 


(Rate gyros in common use today are easily constructed to measure the components of 
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the inertial angular velocity of a vehicle that lie along the vehicle's body axes.) It follows 
that the instantaneous angular velocity of the vehicle can be related to the instantaneous 
rate of change of angular orientation only after considerations of the intermediate 
transformations used. In other words, if @ is defmed as the angular attitude vector, 
a=[¢,4]', and the inertial components of the vehicle angular rate lying along the body 
axes are defined as w= [p,gq, r]’, then & = f(a,q@). 

The details of the nonlinear functional relations involved are provided by viewing 
the rate of change of the rotation y, as a vector quantity lying along the original Z-axis. 
The rate of change of the angle @, 1s viewed as a vector quantity lying along the Y-axis of 
the first intermediate frame, and the rate of change of the angle ¢, is viewed as a vector 
lying along the X-axis of the final (body fixed) frame. Each of these component rates of 
change of angular position has component parts that project onto the final frame and it is 
the sum total of all the components that give the total angular velocity as seen in the final 
frame of reference. Using the required transformations for the rate components from 


each Euler angle, 


p 0 0 d d 
q|=>T@T (OT) 0 |+ (TPT) 6 |+T(9)| 0|=7,4,6.y) 4 |. (2.9) 
r y 0 0 WV 


the body rates are obtained with the result, 


p -ysind+¢ 1 0 -sinod 1¢ 
q\=| wsin@+Ocosd |= 0 cos@ sin 8 6 |.(2.10) 
r W/cosAcos¢ — Osin d 0 —sng  cos@cosd| wy 


Inverting Equation 2.10 yields a solution for the rates of change of the Euler angles in 


terms of the body fixed components of the angular velocity vector, 


d p+qsingtan@+rcos¢ tand ] singtan@ cos@tand | p 

g\= qgcosg@—rsin d = 0 cos@ —sin @ q |. 

yW (qgsin@+rcos¢)/cosé 0 —sin @/cos@ cos¢@/cos@ | r 
(205) 


Notice that for small angular rotations, as expected, 
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d= p; 

0 = 4; 

War, 
as expected. As a note, it should be pointed out that unlike the transformation matrix T;, 
T>2 is not orthogonal, therefore, T.- 51 ¢T,!. 

At this point the kinematic relationships between velocities, as seen in the body 
fixed frame, and the rates of change of global positions and Euler angles have been 
defined. The resulting set of differential equations forms a consistent set in that given a 
set of vehicle velocity data versus time, its position and attitude may be computed. 
Expending and combining equations 2.8 and 2.11, this set of differential equations is 
summarized here; 

X ucos@siny + v(—cos dsiny +singsin @cosy) + w(sindsiny + cos@sin @cosy 
Y ucos@siny + v(cos¢cosy + singsin Osiny) + w(—-singcosy + cos dsindsiny) 


LA —usin@+vsingcos@ + weosgcosé 
d . p+qsingtan@ +rcos ¢tand 
6 gcos¢@—rsing 
Wy (qsing + rcos¢)/cos@ | 


In vector-matrix form, this set of equations can be expressed as; 
aie 0 
z= ax 
0 T, 


Cc: NEWTON’S SECOND LAW 


Since the motion of the vehicle is composed of both translational and rotational 
components, it is necessary to retain the distinction between the vehicle or body fixed 
frame and the global or inertial frame. This is particularly important because the 
dominant forces that act on a submerged body depend on the local motion of the vehicle 
relative to the fluid particles, not on the global velocities. Returning to Figure 2.1, the 


global position, velocity and acceleration vectors of the vehicle’s center of gravity is 


dr. _.d? 
defined as r, =ry + Pg»—& and —2 


7 re respectively. Since the center of mass lies in a 
[ [ 
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body that is both translating and rotating, the total time derivative of rg comes from two 
parts. The first part is from the time derivative as if the body fixed frame was not 
rotating, while the second part is due to the rotation of the body fixed frame. A detailed 
explanation of this can be found in [Greenwood 1988]. 

For a general vector r, in a frame rotating at an angular velocity @ the total 


derivative 1s given as 
an =r+@xr (2.12) 
dt 


where the first right hand side term is a rate of growth part and the second right hand side 
term is a rate of transport part. Using this logic, the velocity of the vehicle’s center of 


mass is given by 


Be 
Gp ato OX Pg =%o + OX Po. (2.13) 


The inertial velocity of the origin of the body fixed frame is expressed as vg: and may be 
written in either global or body fixed coordinates, therefore, 


dX dY dZ 
v,.=Fr,. =[—1+— J + — K ] =[ui+vji+ wk). 2.14 
0 Fo = dt i ee ie 


1. Translational Equations of Motion 


With the coordinate frames defined, Newton's second law of motion, 
fae == (mv >)» May be used to formulate a translational model of the system. The 
f 


global acceleration of the center of mass is derived by differentiating the velocity vector, 


r,, realizing that the center of mass lies in a rotating reference frame. Considering the 


total differential, the global acceleration of the center of mass becomes, 
Tr, =Vo.+ OX Po +OXOX Po +OXV>o.. (25s) 
The translational equation of motion is found by equating the product of this 
acceleration and the vehicle mass, to the net sum of all forces acting on the vehicle in 
three translational degrees of freedom (X, Y, Z). One important factor to recognize is that 
the equation of motion derived in this manner is a vector equation with the components 


expressed in the body fixed frame with unit vectors i, j andk. As discussed earlier, this 
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has been deliberately done because the dominant forces acting on a submerged body in 
motion are developed in relation to the shape of the vehicle and are more conveniently 
expressed in relation to the body axes. Applying Newton's second law results in the 
vector equation, 


feans = NMVo tOX Po tOxXOxXp,t+OXvo}. (2.16) 


Trans 


Ze Rotational Equations of Motion 


To develop the rotational equations of motion, the sum of applied moments about 
the vehicle's center of mass is equated to the rate of change of angular momentum of the 
vehicle about its center of mass. In the practical case of marine vehicles, however, the 
statement just made is modified slightly because it is much more difficult to assess the 
vehicle's mass moments of inertia about its center of gravity (CG), as the CG changes 
with loading. It becomes simpler to evaluate the mass moments of inertia about the body 
fixed frame which tends to lie along the vehicle's axes of symmetry. The inertia tensor 


required to be computed is, 


Lei 
Ds: = eed elite (2.17) 
Pele ibe 
where, 
N N 
l= es dm.(y* Za) and J, =1,, = -> dm, (xy), for example. 
i=] 
The angular momentum of the body is given by, 
h, =1,@. (2.18) 


The total applied rotational moments about the vehicle's reference frame origin is given 


by, . 
| m,, =h, + Po X(mvc). (2.19) 


Differentiating Equation 2.18, the rate of change of angular momentum is found to be, 
h,, =1,@+@xh,, (2.20) 


and the acceleration of the global position vector rq:, is given by, 
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r, =v, +@Xv,.. (2.21) 
Substituting equations 2.20 and 2.21 into Equation 2.19, the rotational equation of motion 
in vector form is given by, 
m,,, =1,@+@x(I,@)+m{p, Xv,.+ Po XOXv,}. (2.22) 
(Again, for a more detailed discussion, the reader is referred to [Greenwood 1988}). 


3 Equations of Motion in Free Space 


At this point, there are three translational equations obtained from Equation 2.16, 
three rotational equations obtained from Equation 2.22, and six unknown velocities 


(v and @). This set of equation written in long form is; 


mlu — vr +wq—X%g(q° +r°)+ yo(pq—F)t+2g(pr+gl=X,, (2.23) 
m[v+ur—wp+Xxo(pgtr)—yo(p’ +r’ )+2Z,6(gr- p)]=Y ‘ (2.24) 
m[w —uq+ vp +x¢(pr-4q)+yo(qr+ P)-Z6(p’ +q I=Z,, (2.25) 
I,p+U,-1,)gr+1,,(pr-q)-1,,(q° -r°)-1,(pqtr)+ oan 
ml yg(w-ugt+vp)—Z,(vtur—wp)|=K, . 
1,g+(U, -1,)pr-1,(qr+p)+1,.(pq-n+1,.(p -r°?)- 

. . (225 
m[xg(w-ug+vp)—Z,6(u-vr+wq@)j=M , 

and 

17+, —-1,)pq-1,(p* -—4°)-1,,(prt+q@)+1,.(qr- p) + 
prot (ie ) Qari Pe) ae mCP Taiatan) (q De (2.28) 


m[xg(v +ur—wp)- yg(u-vr+w@)]=N, 

These preceding six equations are the most general form of Newton's laws for 
rigid body motion used today. They consist of a constant mass and inertia tensor, and are 
formulated in the body fixed frame. This set of equations can be simplified depending on 
the location of the local axes. If the axes are selected to coincide with the vehicle's 
principal axis of inertia, then the terms including the product of inertia become zero. 
This simplification is possible only if the xy-, xz-, and yz-planes are planes of symmetry. 
This is typically not the case. (Practically, most designs are only symmetric about the xz- 


plane. Further simplifications can be made, including locating the origin of the body 


20 


fixed frame at the vehicle's center of gravity, however, these assumptions are not 
physically realizable.) 

The previous nonlinear system, equations 2.23-2.28, can be written in a compact 
vector form. Using the vector x the state vector of body fixed velocities, and defining 


Mpg as arigid body mass matrix including translational and rotational] inertial elements, 


m 0 0 Oi get) = 
O m OF iz ae LX 
0 OR leery ee eee 
M ee = OM ize eae” los Ti | 
Yi ee OX eed 3 
= OVE ee es 


and Crp(x) as a State dependent matrix containing the rigid body coriolis and centrifugal 


terms, 
Crp (x)= 
0 =i) mq mM Y6q + Z¢r) —MX¢q MX er 
mr 0 —mp TGP m(Z¢rt X¢P) —Mer 
-mq mp 0 ~ Mg Pp — mq m(X¢Pt+Yeq) |; 
—M(Yoq+ Zor) MYcP MZ¢p 0 alg) Opti lr ce pie 
MX oq —m(Zor+XxGP) MZcq hd td ap ar 0 =o leGa) ep 
mXor myor SING Ped) ete Pig. nel gp 0 


the equations of motion can be written in vector form as, 
M ~x+C (x= fj? (2.29) 
The right hand term in Equation 2.29 is the vector of external forces and moments 
outlined in Table 2.1. These forces and moments come from hydrostatic and 
hydrodynamic forces due to gravity, radiation and excitation, viscous damping (lift and 
drag), and control inputs. The omgin of these forces and thew application to the 


developed equations of motion will be discussed in detail in following sections. 


D. RESTORING FORCES AND MOMENTS 


In hydrodynamic terminology, the gravitational and buoyant forces are called 
restoring forces. The weight, W, and buoyant, B, forces that act at the centers of gravity 
and buoyancy must be defined from static analyses. For submerged bodies the weight 


and buoyancy force vectors do not change with vehicle attitude. Assuming that weight 
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and buoyancy are fixed in relation to the body fixed frame, then the gravitational and 


buoyant forces may be expressed as, f,, =O1+0OJ+WK, and f, =O1+0J-BK. 


Since the weight and buoyancy terms in the applied forces act in the global vertical 
direction, they must be transformed into components in the vehicle fixed frame before 
they can be added into the equations of motion. Returning to Equation 2.4, it can be seen 
that the components acting along the vertical vehicle fixed frame are the third column of 

the transformation matrix. Therefore, the net vertical force components become, 

— sin 8 
f,, =(W — B) cos@sin ¢ |. (2.30) 
cos@ cos @ 

The weight portion of the vertical force acts at the center of gravity of the vehicle, 
while the buoyancy portion of the vertical force acts at the center of buoyancy. Because 
these forces act in locations away from the body center they create a moment about the 


body center given by, 


—sin@ — sin 8 
m, =Wp, X| cos@sin@ |— Bp, X| cos@sing |. (Za 
cosOcos@ cos@cos@ 


This moment will be non-zero even if W and B are equal, since Pg and pp are 
usually not collocated. For static stability it 1s necessary to locate the center of gravity 


below the center of buoyancy. The total vertical force vector can be written as, 


—(W-B)sin@ 
(W —B)cos@sin @ 
: ue 7 (W —B)cos@cos@ =a 
fo|3" |- (y,W —y,B)cos@cos@—(z,.W —z,B)cos@sin F(). aie 


—(x,W -—x,B)cos@cos@—(z,W -z, B)sind 
| (xgW-x,B)cosOsin 9+ (y,W-y,B)sin @ 


and is added to the right hand side of the equations of motion. 


Ze 


E. WAVE THEORY AND HYDRODYNAMICS 
ie Linear Wave Theory 


The simplest free surface wave formation, which nevertheless has great practical 
significance, is the plane progressive wave system where the water column is modeled as 
an inviscid, irrotational fluid in a gravity field. This motion is two dimensional, (x, 2), 
sinusoidal in time with angular frequency @ and propagates with a phase velocity c, such 
that to an observer moving with this speed the wave appears to be stationary. A 
Cartesian coordinate system is adopted, see Figure 2.2, with z = 0, the plane of the 
undisturbed free surface (still water level) and the z-axis positive upwards. The vertical 
elevation of any point on the free surface may be defined by a function z = 7(x,t). With 
these requirements, the free surface elevation must be of the general form 

7(x,t) = Acos(kx — at), (2733) 
where the positive x-axis is chosen to coincide with the direction of wave propagation. 


Here, A is the wave amplitude, and the parameter 


en (2.34) 


Cp 


is the wavenumber, the number of waves per unit distance along the x-axis. Clearly, the 
wavenumber can also be written as, 


k still (2.35) 


L 
where the wavelength L, is the distance between successive points on the wave with the 
same phase, see Figure 2.2. The solution of this problem is expressed in terms of a two 
dimensional velocity potential, which must satisfy Laplace's equation 
V*o=0, (2.36) 
and appropriate boundary conditions. Furthermore, the velocity potential, ¢, must yield 


the wave elevation given by Equation 2.33 from 


Vea ae GL LHe (2.37) 


Zo 


a Progressive Wave 


L Direction of motion 
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Figure 2.2, Monochromatic Progressive Surface Gravity Wave [Rahman 1994] 


Equation 2.37 is the linearized dynamic boundary condition on the free surface and is an 
expression of the fact, through Bernoulli's equation, that the pressure on the free surface 
must be the same as the ambient atmospheric pressure. An appropriate boundary 


condition on the sea bottom is 


eae (2.38) 
Oz 


i.e., the bottom at depth H is a rigid impermeable plane. Finally, the free surface 
boundary condition is 
d°g a6 
—+g—-=0, at z=0. Zoe 
ar ° a = 
Equation 2.39 is a combined dynamic and kinematic surface boundary condition. 
The dynamic condition was discussed earlier, while the kinematic condition simply 


states, 


dn oO 

see. (2.40) 
that the vertical velocities of the free surface and fluid particles are the same. Combining 
equations 2.37 and 2.40, Equation 2.39 is obtained, ignoring the small departures of the 
free surface 7) from the horizontal orientation of z = 0. 


From the requirements of the problem, it is clear that the velocity potential @ must 


be sinusoidal in the same sense as Equation 2.30; therefore a solution of the form 
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Q(x, z,t) =sin(kx — at) F(z), (2.41) 
is attempted. Substituting Equation 2.41 into Equation 2.36, F(z) must satisfy the 
ordinary differential equation 

a? F(2) 
dt? 
throughout the domain of the fluid. The solution to Equation 2.42 satisfying the bottom 


—k*F(z)=0, (2.42) 


boundary condition is 

F(z)=Acosh(k(z+f)). (2.43) 
Substitution of equations 2.41 and 2.43 into the surface boundary condition, Equation 
2.37, yields an important relationship between the wavenumber k and the frequency @ 

a = gk tanh(kH), (2.44) 
which is called the dispersion relationship. The surface elevation 7 follows from 
Equation 2.37 as, 

7(x,t) =acos(kx — at), (2.45) 


with the amplitude a, given by 


SET (2.46) 
g 


Substitution of equations 2.43 and 2.46 into the velocity potential function, Equation 
2.41, yields 


ag cosh(k(z + #7)) 


in( kx — ar) . 4 
@  cosh(kH) a se ee 


Q(x, Z,t) = 


An underlying assumption associated with potential flow theory and Laplace's 
Equation, is that the velocity field can be expressed as the gradient of a velocity potential 


function. This allows the expressions for the velocities in the x and z direction to be 


given as 
w= 
im 
; (2.48) 
a 
dz 


Using Equation 2.48, and the linearized form of the Bernoulli's equation, 
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dg 
—— 2.4 
p Po pPgz (2.49) 


the expressions for the fluid velocity and pressure fields are 


png Bk cosh(K(z + H)) 


cos(kx — at) 
@  cosh(kH) 
ate eee ee | (2.50) 
@  cosh(kH) 
_ cosh(k(z + A)) 7 a 
p= Pea snk) cos(kx — at) — pgz 


It can be seen from Equation 2.50 that the trajectories of the fluid particles are elliptical. 

There are several simplifications that may be made to the above-derived 
expressions for the cases of shallow (long waves) and deep (short waves) water. The 
shallow and deep water ranges correspond to H/L < 7/10 and AVL > f, respectively. 
Over these ranges approximate expressions may be substituted for the hyperbolic 
functions that have been encountered. Table 2.2 summarizes these results. Figure 2.3 
depicts the comparison of water particle velocity between long and short waves. The key 
points that should be observed are that shallow water waves are non-dispersive, i.e., the 
vertical component of the wave particle velocity is linear in depth, and that the 
classification of shallow water depends on the ratio of water depth to wavelength. 

As an example, consider a one-meter high, monochromatic wave with a ten- 
second period (0.1 Hz), propagating in water six meters deep. The deep water and 
shallow water wavelengths are 156 meters and 76.7 meters, respectively. The ratios of 
water depth to wavelength are 0.038, for the deep water case, and 0.078, for the shallow 
water case. Referring to Table 2.2, it can be seen that neither the deep water nor shallow 
water simplifications may be used. In fact, the water depth must be reduced to 
2.45 meters before the shallow water equations are valid. Conversely, if the water depth 
remains at 6 meters, then the wave period must become greater than 15.6 seconds or the 
wave frequency less than 0.064 Hz. This indicates that only low frequency waves 
typically qualify as shallow water (non-dispersive) waves which becomes an important 


consideration in the proper modeling of the wave induced disturbances. 
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Table 2.2, Summary Of Small Amplitude Linear Wave Equations 


As a Side note, it should be pointed out that the wave height has no bearing in 
determining whether a wave is classified as long or short. The wave height is significant 
in order to determine the subsurface water particle velocities, and when the wave breaks. 

The plane progressive wave described so far is a single, discrete wave system, 
with a prescribed monochromatic component of frequency @ and wavenumber k, moving 


in the positive x-direction. More general wave motions, which are not monochromatic, 
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Figure 2.3 Particle Orbits And Variation Of Particle Velocity Amplitude With Depth 
[Sarpkaya 1981] 


can be obtained by superimposing plane waves of different frequencies and 
wavenumbers. The elevation of the sea surface 7(t) can thus be described as the 


superposition of an infinite number of sinusoids of the form: 


n(t) = Sa, cos(k,x-@,t+@,) =>n, ; (2:5) 
n=l n=l 
The equations for the horizontal and vertical velocities, and the dynamic pressure are, 
WOE 5] SED, (2.52) 
w(t)= > et, ; (2353) 
ptt) = pa 3) rsoshta +2, — PRZ; (2.54) 
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respectively, where the parameters/variables in equations 2.51-2.54 are identical to those 
in the monochromatic case. 


Zs Hydrodynamics 


[Newman 1977] has shown that the total velocity potential ¢, may be written as 


~=Oy, +9, +,, the linear sum of three components. These three components come 


from wave, diffraction and radiation potentials where, 

e gy is the incident regular wave velocity potential; 

© @p is the diffraction potential caused by reflection when the vehicle is considered 
restrained from motion; and 

e rp is the radiation potentials in 6 DOF caused by forcing the vehicle to oscillate with 
wave excitation frequency, when there are no incident waves present. 

Linear wave theory implies that the wave induced forces and moments acting on a 

vehicle come from the superposition of the radiation induced forces and moments and the 

excitation forces and moments. 

Radiation induced forces and moments act on the vehicle when the vehicle is 
forced to oscillate with the wave excitation frequency without incident waves present. 
The forces due to wave radiation are classified as restoring, added mass and potential 
damping forces. 

Excitation forces and moments act on the vehicle when the vehicle is restrained 
from motion and incident waves are present. These forces and moments caused by wave 
excitation are classified as Froude-Kriloff (F'K) and diffraction forces. The FK forces and 
moments are found by integrating the pressure distribution over the vehicle cause by the 
undisturbed wave field, while the diffraction forces and moments are determined by the 
pressure distribution created when the waves are reflected from the vehicle. 

As a general summary, the unforced equations of motion for any stationary body 
in waves, can be given as 


(M gp + A(@))7] + B(@) + C7 —-(M x + A(@))7, —B(@)n, —C, =0, (2.55) 


where 77 is a vector of linear and angular displacements. Stationary in this sense implies 


that the vehicle is non-maneuvering and that its only motion is that caused by the body’s 
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interaction with the waves. As seen, the added mass and damping matrices, A(@ and 
B(q@ respectively, are functions of the incident wave frequency. However, if the vehicle 
is moving with some forward velocity U, the wave frequency, @ is not the frequency 
encountered by the vehicle. The encounter frequency, @, can be expressed as 


@,=a+kU cos PB. (2.56) 


In this equation, which is based on the Doppler effect, f is the heading angle between the 
vehicle and the wave propagation direction and k is the wavenumber. More detailed 
information on the Doppler effect maybe found in Appendix B. 

As a simplification, if the body is small in comparison to the wavelength, it is 
totally submerged and neutrally buoyant with ROM sehen mass distribution, then 
Equation 2.55 can be simplified to yield 

(Ms + A(@,))#, + B(@, 7, +C7, =0, (2.57) 


where 77, =7] —7], , 1s the relative velocity of the fluid over the vehicle. It is necessary to 


point out that the above equation was based.on the assumption of inviscid flow. Due to 
this fact, viscous damping terms (skin friction and drag), must be taken into account for a 
vehicle operating in a real fluid, to complete the model. 

In 1950, Morison developed an expression for the horizontal force on a 
cylindrical pile subjected to waves. His work showed that the elemental force dF ona 


vertical strip of a cylinder dz, may be written as 
mD* sea 
dF = p——deC yt, +5 pCpDdzu ; |u|. (2.58) 


In this expression, coined ‘“‘Morison's equation’’, D is the characteristic diameter, u, 


and uy are the horizontal component of the undisturbed fluid acceleration and velocity, 
respectively, and Cy and Cp are mass and drag coefficients which may be determined 
experimentally. This equation has two parts, the first term representing an inertia force 
proportional to the accelerating fluid acting on the pile (a mass force), and the second 
term, a nonlinear drag term proportional to the sign squared fluid velocity (a drag force). 
In practice, Morison's equation can only be applied to small volume bodies. By 


small volume bodies, it is meant that the characteristic cross sectional dimension of the 
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body is small relative to the wavelength [Morison 1950]. For vertical piles, small volume 
applies only if L>5D, where D is the pile diameter, refer to Figure 2.4. As shown, for 
small vehicles operating in shallow water, inertia and drag forces are dominant, where as 
reflection and diffraction effects can be considered unimportant. 

It is known that the ratios of wavelength and wave height to the characteristic 
diameter are key parameters in predicting the load regime of waves acting on a structure 
[Faltinsen 1990]. These regimes are also depicted in Figure 2.4. 

For a stationary object in a simple harmonic flow, the total force can be expressed 
as 

F, =F, sin cx|sin cot| + F, cosa, (2559) 
where Fp and F; represent the maxima of the drag and inertia force components, 
respectively. [Dean 1984] divided the flow force regimes into two areas; one where the 
inertially derived component dominated the total force, and one where both drag and 


inertial effects are important. This expression is represented as 
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Figure 2.4 Load Regimes 
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The significance of the above expression is that the maximum force on the body is not 
affected by additional drag force until the amplitude of the drag is at least one half that of 
the inertia force. For oscillating flows, such as those caused by waves, while even small 
amounts of drag may be important when considering the shape of the load function on a 
stationary body, the peak amplitude of the force is only affected when the drag 
component is greater than one-half the inertia force. Extensive experimental verification 
of Morison’s approach to force modeling, and the evaluation of the frequency dependent 
nature of the drag and added mass coefficients was given by [Sarpkaya 1975]. 

When the wave field causes motion of the water particles at the vehicle's 
operating depth to be of the order of one diameter or less, it is expected that the 
predominant hydrodynamic force on the AUV due to the wave disturbance would be 
inertial in nature. Since the AUV operates below the surface, it is not the wave height to 
vehicle diameter ratio that is of concern, but more appropriately the double amplitude of 
the water particle motion at the vehicle operating depth compared to the vehicle 
characteristic length of interest. While this analogy is approximate in nature, it does 
provide a means to predict which hydrodynamic forces may be of concern when 
estimating the total load on a vehicle. These above concepts may be used to estimate the 
dominant forces acting on an underwater vehicle subject to waves, and therefore assist in 


the sizing of the propulsion system. 


F. HYDRODYNAMIC FORCES AND MOMENTS 


Hydrodynamic forces and moments are the result of body/fluid interactions. The 
forces and moments on the body arise from the modification to the pressure distribution 
summed around the surface area of the body. This modification to the pressure field can 
only arise from relative velocity and acceleration between the body and fluid. Therefore, 
for the purposes of this discussion it is necessary to re-define the body fixed velocity 
vector x, in terms of a relative body fixed velocity vector x,, where 
MV Wee ae ry" Also at this time it is convenient to define a globally based fluid 


velocity vector Uy, where Ur = [U,, Vp, W,, 0, 0, Ol. Since it 1s assumed that the fluid 


a2 


velocity is irrotational, no changes to the angular rate terms are necessary in the body 
fixed velocity vector, and no angular rates are present in the fluid velocity vector. 


1. Radiation Induced Forces and Moments 


a) Added Mass 


Like the rigid body kinematics, it is desirable to separate the added mass 
terms into terms which belong to an added mass matrix M4y and a matrix of Coriolis and 
centrifugal terms C4y(x). For underwater vehicles this implies that the added mass 
forces and moments can be written as: 

Faw =—M ay %, —Cay(®,)x, (2.61) 
where fau=[Xa, Ys,Za,Ka,Ma,Na]'=[fa, ma la is the total added mass force and moment 
vector. 

The added mass terms represent the inertial reaction of fluid particles 
surrounding the submerged body that are accelerated with it. Any motion of the vehicle 
induces a motion in the otherwise stationary fluid. In order to allow the vehicle to pass 
through the fluid, the fluid must move aside and then close behind the vehicle. As a 
consequence, the fluid passage possesses kinetic energy that it would lack if the vehicle 
was not in motion. 

[Lamb 1932] gives the following expression for the fluid kinetic energy, 
E,, which may be expressed in a quadratic form of the body axis velocity vector 


components; 
] 0@,, lf 
ee =-5 P|, os 45 =>, Me (2.61) 


In Equation 2.61, Mam is a 6x6 added mass matrix. For a rigid body moving in an ideal 
fluid the added mass matrix is symmetrical, 1.e., May = Mon In a real fluid these 36 
elements may all be distinct. [Wendel 1956] has shown that the numerical values of 
Wiaee mass in a real fluid are usually in good agreement with those obtained from ideal 
theory. For a body possessing vertical plane symmetry only, and applying Mam = Mam . 


the added mass matrix is written as 
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The notation of [SNAME 1950] is used in this expression. This notation indicates the 
degree of freedom on which the hydrodynamic added mass force acts, as well as the 


cause of the force. As an example, Y,u is a force acting along the body fixed y-axis due 


to an acceleration u in the x-direction, and can be thought of mathematically as 


Y,= = This definition implies that the hydrodynamic derivatives corresponding to the 


“Out 
diagonal of the added mass matrix will all be negative. 

The added mass terms are obtained from potential theory. This theory 
assumes an inviscid fluid, no circulation and that the body is completed submerged in an 
unbounded fluid. The last assumption is violated at the seabed, near underwater objects 
and at the surface. [Milne-Thomson 1968] has shown that the expressions for the force 
and moments resulting from added mass effects can be found by applying Kirchhoff's 


equations, 
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ie rox oe =m 
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in vector form. Expanding Equation 2.63, the expression for the added mass terms 
associated with the x-direction is 
X, =X ,u+X,,(w+uq)+X ,g+Z,wq+Z,q° 
, (2.64) 
=T VW I pa 
(Derivation of the added mass terms associated with the other degrees of freedom is left 


to the reader.) Many of the added mass derivatives contained,in the general form are 


either zero or mutually related to another term when the body has various symmetries. A 
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more detailed discussion of added mass terms applied to an underwater vehicle, is found 
in [Humphreys 1978]. 
Extracting the added mass derivatives corresponding to the velocity 


coupling terms from 2.63 yields; 


> 
© 
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e) 
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where, 
Cis =-X,u—-Z,w-Z,q Cy,=X,u+X,wt+X gq 
Cig =Y¥,v+Y,pt+Y,r GRASSO SISOS IN Se 
Cy =X, utZiwtZ,g Cy =X ,ut+Z,wt+M,¢q. 
Cx, =-XU-X,w-X,q Co =-Y,v-K,p—N,r 
Cy =—-Y,v-Y,p—Y,r 


b) Wave Radiation or Potential Damping 


The contribution from the potential damping terms compared to other 
dissipative terms like viscous damping terms are usually negligible for underwater 
vehicles operating at great depth. Nevertheless, underwater vehicles operating in shallow 
water near to the free surface should consider potential damping effects, especially, those 
underwater vehicles that tend to have a non-streamlined body, 1.e., vehicles build with 
sensor and equipment mounted in such a manner that the equipment causes the vehicle 
not be streamlined. The linear potential damping can be modeled as 

f, =—-D,x,. (2.65) 
The linear damping matrix Dg is a positive definite matrix of linear damping coefficients. 
These linear damping terms are small when compared to the viscous forces, and therefore 


are often included in the viscous drag forces. 
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Da Excitation Forces 


When applying potential theory, the fluid motion was assumed to be urotational. 
This implies that only the linear velocity components, Uy, of the fluid are considered 
when determining the excitation forces. In linear theory, the wave induced forces and 
moments acting on a vehicle are considered to be the sum of the radiation induced forces 
and moments and the excitation forces and moments. For nonlinear theory, this is not the 
case. The forces and moments due to radiation and diffraction are nonlinear functions of 


the relative velocity and acceleration between the vehicle and the fluid. 


a) Froude-Kriloff Forces 
The Froude-Kriloff force and moment vector can be expressed as 
Sirk = M px, , (2.66) 
where Mrx can be interpreted as the Froude-Kriloff (FAK) mass matrix and u y 3S the fluid 


acceleration expressed in body fixed coordinates. Coriolis and centrifugal terms will not 
appear in the general expression for the FK forces and moments since it was assumed that 
the rotational fluid motion was zero. This mass matrix can be determined by computing 
the mass of the fluid displaced by the vehicle and substituting this value into the rigid 
body mass matrix in place of the vehicle mass. Also, it should be recognized that this 
excitation force acts at the center of buoyancy not the center of gravity. The mass and the 


moments and products of inertia for the PK mass matrix are 
N 
m = pV =B/g, T,, = >._,.dm,(y? +2”) and I, =1,, =-> dm, (2x9). 
i=] 


Substituting these values into the rigid body mass matrix yields the FK mass matrix, for a 


small volume completely submerged body. 


m 0 0 0 iz ee TS 
0 m 0 =z 0 MX p 
0 0 m iy pe 0 
M rx = 0 —Nlz pany a ly |e 
MZ, 0 tot oe oe l,, 
Pity, One, 0 a ih ik | 


[Sarpkaya 1981] and [Newman 1977] have shown that the body fixed 
inertia force to which a symmetric body moving in an unsteady flow field is subject, may 
be written as 


. ou , Ou 
Te ON It) eee (Uy) |= ON Ci. (2.67) 
ot Ox 


In this equation, it is important to note the presence of a buoyancy-like force that is 
proportional to the displaced volume of the body as well as the presence of some 
convective terms . The terms that include the displaced volume, when grouped together, 
represent the added mass and Froude-Kriloff forces. The convective terms represent the 
forces from the spacial changes of the unsteady flow field over the body, and for large 
bodies, may be significant. [Silvestre 1998b] has shown that by comparing the force 
contribution of the convective terms to the total inertia force exerted on a rigid body 
subject to wave disturbances, that the forces due to the convective terms are small and for 
all practical purposes may be neglected. 


3. Viscous Damping Forces and Moments 


a) Drag Force 


The drag effects on an underwater vehicles are mainly caused by 


2 Linear skin friction due to laminar boundary layers; 
3 Quadratic skin friction due to turbulent boundary layers; and 
3 Quadratic drag due to vortex shedding (Morison's equation). 


The viscous damping forces and moments will be functions of the relative fluid motion. 
In the range of Reynolds numbers in which underwater vehicles typically operate, flow is 
turbulent, therefore the drag force is approximated by the square law resistance arising 
from Morison's equation. Referring to Equation 2.55, the quadratic drag force in the 


x-direction can be expressed as 








f= ; pC, Au, (2.68) 


u, 


where A is the projected cross-sectional area, Cp is the drag-coefficient based on the 


representative area, p is the fluid density and u, is the relative longitudinal velocity. 
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A generalization of Morison's equation could be to use a truncated second 
order Taylor series expansion to describe the viscous damping in 6 DOF. This suggests 
that the viscous damping forces could be written as 


f, =-D,x, —D,(x,)—higer order terms . (2.69) 
This expression is a combination of linear and quadratic drag components. The linear 


portion D,x,, can be modeled as a positive definite matrix of linear damping coefficients 
multiplied by the body fixed relative velocity vector, while the quadratic portion D, (x, ), 


is a nonlinear matrix incorporating the contributions due to skin drag and vortex 
shedding. 

It is quite complicated to determine the hydrodynamic coefficients 
associated with the quadratic portion, especially for high angles of attack. However, a 
simple, but fairly accurate method of modeling these viscous drag forces is to include the 


linear drag components associated with the diagonal of D,x,, and to use a cross flow 
integral to account for the D,(x,) terms. The general cross flow drag expression for a 


body of revolution in the y-direction (sway) is given by 


i? 
=-5P [ (Cp, Dn + xr)? + Cy, Dw - xq) 2a. (2.70) 


Vv 
ae U 5 (x) 


drag 


In this expression, the minus sign is present because the drag force opposes the motion, 
and the cross flow velocity is given by 

Uy (x) =[(v + xr)? + (w—x9)7 J”. Qa 

Equation 2.70 is valid for a vehicle that has a hull form consistent with a 

body of revolution, however, for underwater vehicle's with different shapes this equation 

must be modified. As an example, consider the box shape hull form of the NPS Phoenix 

AUV [Marco 1996]. The cross flow drag expression in the sway direction for this shape 


vehicle becomes 


Y 


drag 


1/2 
= -=p es A(x)(vt xr)(v + xr) . (27) 


-l/2 
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As depicted in Figure 2.5, this expression reflects the drag force over the entire range of 


angles of attack. 

b) Lift Forces 

Lift is a force acting on a body in a direction perpendicular to that of the 
flow of the fluid. As relative motion is created between the underwater vehicle and its 
fluid environment, the vehicle body experiences lift forces similar to those experienced 
by an airfoil. The expression for the lift generated by an airfoil is given as 


dC 
poy = an (2.73) 
p da \u, 


where A is the projected cross-sectional area, Cz is the lift-coefficient, @ is the local angle 
of attack, p is the fluid density, and wu, and v, are the relative velocity components in the 
body fixed x and y directions. See [Hoerner 1975] for further detail. 
The lift forces can be modeled as a state dependent matrix multiplied by 
- the respective relative velocity and is given by 
f, =—D,(x,)x,. (2.74) 
The parameters/coefficients used in this matrix, like the drag coefficients are difficult to 
predict and vary with the shape of the body and the angle of attack of the fluid. As with 
the drag forces, determination of these constant coefficients is typically done using a 3-D 


CFD solver, with the constant coefficients validated over a range of angles of attack. 


G. CONTROL FORCES 


Small underwater vehicles are usually maneuvered with thrusters and control 
surfaces. Thrusters are effective only at low forward vehicle speed due to the fact that 
the action of the thrust is a nonlinear function of the relative velocity of the vehicle. 
However, control surfaces are effectively used in maintaining heading as well as trim and 
depth changing maneuvers. The reason for this is that the force generating capacity of 
the control surfaces is dependent on the speed of the vehicle (the lift force is proportional 


to the square of the velocity). The control forces and moments can be described by 


Sc oa B(x, > UL control | es (2775) 
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Figure 2.5 Normalized Drag Force vs. Angle of Attack 


where Ucontrol 1S aN input vector and B(Xpy Ucontroi) 1S a State dependent input matrix. This 
input matrix contains the necessary coefficients to model the forces developed by the 
control actuators. The total force/moment vector fc, is the sum of the propulsion force 
vector, fp, due to thrusters and propellers, and the actuator force vector, fs, due to fin and 
rudder deflection. 


1. Propulsion Forces 


The derivation of a steady-state hydrodynamic model for propellers operating in 
an incompressible fluid can be found in many introductory fluid texts see [Lewis 1988, 
Newman 1977 and White 1986] for examples. These models are based on large open 
propellers. The use of small thrusters for control of underwater vehicles is an area of 
current research. This is because the vehicles are small, require fast response and are 


required to conduct dynamic positioning maneuvers. [Yoeger 1991] developed a lumped 
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parameter model that improved on the popular notion that for a given unit the thrust and 
input torque are related to the square of the propeller rotational rate and the angle of 
advance. He also showed that by accounting for thruster dynamics, improvements in 
position control could be obtained. Yoeger's work introduced the concept that fluid 
momentum considerations resulted in a time lag in the thrust response to a step input. 
Although his work improved the modeling of thrusters, it left room for significant 
improvements. 
| [Healey et al., 1995], improved on Yoeger's work by providing a generic thruster 
- model that considered propeller thrust and torque as a mapping linked to lift/drag force 
variations caused by changes in the local angle of attack of the propeller blade. They also 
were able to associate the lags and overshoots in the thrust response to lags in the 
development of the local angle of attack and dynamic development of the blade pressure 
distributions. Research by [Whitcomb 1995] and [Bachmayer 1998] has provided further 
experimental validation that the four-quadrant model proposed by Healey is valid for 
small ducted thrusters. 

The model proposed by the researchers in the above paragraph can be best 


expressed as a first order differential equation of the form 


+ F rl : (2.76) 








poe Re n 
t t 


where F is the propulsive force imparted to the vehicle. The parameter 7 is the time 
constant associated with the force lag, the second term containing ¥, is a thrust reduction 
cause by the change in local angle of attack as the propeller advances through the water 
and the last term containing f, is the thrust to rotation rate mapping. In standard 
propeller design terms, the parameter ycan be related to J, while {can be related to Kr. 

Referring to Figure 2.6, it can be seen that the thrust coefficient Ky is a function of 
the propeller speed of advance J. The value of the non-dimensional thrust coefficient at a 
zero speed of advance is associated with the bollard pull condition. This relates the thrust 
to the rotational rate of the propeller as 


T =K,pD*nln\, (2.77) 
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where D is the diameter of the propeller and n is the rotational rate. 
For any other operating condition, it can be seen that the thrust coefficient is 
reduced as the propeller moves through the fluid with some speed of advance. The thrust 
for this operating condition is given by 
T =K,pD*n|n|+y,JpD* nln], (2.78) 
where % is the slope of the Ky curve at the particular operating condition of interest 
(negative over the entire curve), and J is the non-dimensional speed of advance. 
Substituting the expression for J, shown in Figure 2.6, the parameters f and y from 
Equation 2.76 can be related to the parameters in Equation 2.78 as 
B=K,pD° 


- (2.79) 
Y=Y,pD 
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Figure 2.6 Thrust and Torque Coefficients versus Angle of Attack [Lewis 1988] 


The force/moment vector caused by the propulsion forces due to the set of 


thrusters and main propellers for the NPS Phoenix AUV is, 
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where the x and y distances are all positive values measured from their respective plane 
of origin. 


Ze Actuator Forces 


Control of a small AUV at anything other than slow speed must be accomplished 
by using control surfaces since the effect of non-propulsion thrusters decreases as the 
forward speed of the vehicle increases. These control surfaces are comprised of 
fins/planes and rudders. The forces and moments exerted by these actuators are derived 
from airfoil theory and are composed of lift and drag components. With the exception of 
the longitudinal direction, the force/moment applied to the vehicle is directly proportional 
to the amount of angular deflection of the control surface. In the longitudinal direction 
the force exerted by a deflection of a surface amounts to an additional drag force on the 


vehicle. The vector of forces and moments caused by these surfaces is 
[(X 45, Om +X q5, Osp UI+(X 5, Op +X 5,05, Ur + X,3 OUV+... 
(X55, +X 5,5, uwt (X55 (5; +65,)+X55,5, +X, 5 Om ulull 
(Y; 0, +Y5 ,, yulu| 
f= (Zs, 54 +Z5, Sy, )ulu| 231) 
0 
(M Os O5 i Ms, Op ulu| 
(Ns 6,,+Ns, 6,, )ulul 
The modeling represented in Equation 2.81 is for a vehicle equipped with a 
standard fin arrangement. By this it is meant that each pair of control surfaces, stern 
planes, bow planes or rudders, move together and are not independently controlled. If the 


vehicle is equipped with either an X-brace configuration[Humphreys 1994], or the 


actuators in a standard configuration are allowed to be moved independently, then the 
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modeling of these control forces will change. Specifically, there will be the ability to 
actively control roll and there will be some level of redundancy for each of the other 
control modes. For a detailed description of the origin of these control forces, the reader 


is referred to [Hoerner 1965, Hoerner 1975 and Lewis 1988]. 


H. 6DOF EQUATIONS OF MOTION 
1. No Fluid Motion 


Using Equation 2.29, and substituting the external forces discussed in the 


previous sections, the DOF EOM can be written as 
M ggX +Cgg(X)X = SG + fam St ny reat dy elfen oe (2.82) 


Expanding Equation 2.82 by substituting the appropriate expressions for the 
force/moment vectors, and recognizing that without fluid motion x = xr and us = 0, the 
following equation is obtained 

iM ppt tO p(X) oath (2) 4M Cee 


(2.83) 
D,X+D,x+D(x)+ D,(X)xX = BUX, U contro: 


control 

As seen, Equation 2.83 contains a mixture of coordinates; both body fixed and 
global. To use a system of equations expressed in this form for simulation studies or 
control system development, this system must be augmented with the relationships 
between the various coordinate frames. The link between the global and body fixed 
coordinates is accomplished by augmenting Equation 2.83 with equations 2.8 and 2.11. 
This system is represented by 


POE Sek, pe aslC es iar(  pehie ar 


[D, + D, + D,(x))x + Di (X) + F(Z) = BUX contro: contra (2-84) 
Z = g(x,Z) 
Ze Modifications To Account For Fluid Motion 
In the case where fluid motion is present, Equation 2.83 is represented as 
M ppXt+C,p (x)X+F(Z)+M,,,x,+C, (x,)x,+D,x, (2.85) 


—-M xu, + Dx, +D,(x,)+D,(x,)x, = BCX, .U consrot U control 


This equation, as was Equation 2.83, is a mixture of various coordinate frame 
variables; body fixed, body fixed relative, global and fluid. To solve this system of 
equations, the equations must be expressed in variables that can be related to each other. 
Since, for this case, the vehicle is considered to be in an unsteady fluid referenced to the 
global frame (refer to Figure 2.1), the logical choice. for variables is body fixed relative 
xr, and global, z. 

Remembering that relative velocity is defined as xr = x — u,, and that the body 
fixed fluid velocity can be expressed as 


T, 0 
“=| lu, =10,. (2.86) 
0 T, 


Equation 2.85 can be modified so that it is expressed in body fixed relative and global 
variables. Manipulating Equation 2.85 results in 
[M pe +M ay )x, ms [Cra(X, ) TCT Ae IX; 7 [D, 5 5 D, (x, Date D, ee Ix, sez) 
d(TU , ) d(TU ,) 


=-M,, cee TU, +M rx Tp Br Heoni ju 


control 

(2.87) 
Looking at Equation 2.87 and recalling that the fluid was defined as irrotational, it can be 
seen that the Crg(us) term on the right hand side of the equation is zero. The other two 
terms, on the right hand side, containing the rigid body mass matrix Meg, and the 
Froude-Kriloff mass matrix Mrx, are excitation forces resulting from the fluid motion. 

As seen in Equation 2.87, the Froude-Kriloff excitation forces and moments are 
functions of the weight and buoyancy mismatch (W-B), the fluid velocities and fluid 
accelerations, expressed in body fixed values. For a neutrally buoyant vehicle, where 
W=B, the wave excitation forces do not present themselves in the translational equations 
of motion, however, they still provide excitation moments to the rotational equations. 
The reason behind this is because the fluid components (acceleration and velocity) act at 
the vehicle ene of buoyancy, while the body inertial acceleration reaction force acts at 
the vehicle’s center of mass. 

At this point the 6DOF EOM representing the vehicle dynamics has been 


modified to account for a moving fluid by representing the body fixed velocities in 
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relative terms, Equation 2.87. However, as with the system of equations for the case of 
no fluid motion, the system in Equation 2.87 must be augmented to provide the necessary 
link between the global and body fixed velocities. In order to account for the fluid 
motion, either wave induced or steady current, Equation 2.8 is modified and is 


represented as 


X u, U, 
Y |=7,"(9,6,w)| v, |+1V;, |, (2.88) 
Z Ww, Ww, 


where the fluid velocity is represented in the global frame. Combining equations 2.87, 
2.88 and 2.11 into a system of equations results in the necessary equations to describe the 
motion of a small underwater vehicle subject to shallow water waves. 

It was shown earlier in this chapter that the total time derivative is composed of a 
rate of growth term and a rate of transport term. In the case of the time derivative of the 
rotation transformation matrix J, since there is no translation, the time derivative can be 
expressed solely as @T. [Healey 1992a] and [Fossen 1994] have shown that the 
cross-product operation (@x<), can be represented as a matrix by defining the screw 


symmetric matrix, 


O -r gq 
S(@)=| r 0 —pil. (2.89) 
Ga? 


The elements of S(@) are based on the consideration that both vector-matrix quantities in 
the cross-product operation are represented in the same coordinate frame. In the case 
where J is the transformation matrix from global to body fixed coordinates, the 
expression for the rotation matrix time derivative is 
T= : ‘(@r, 0 | (2.90) 
0 L 
S'(q@) is the transpose of the matrix represented in Equation 2.89 since T is not in the 
same coordinate frame as @ Using Equation 2.90, and the fact that the coriolis matrix, 


Cra(uy), is null, the vehicle dynamics equations expressed in matrix form are 
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PY ety aC nn ( Careers) xu (Dex) + Dx ks 
=(M x — M gp [TU , +TU y+ BOX, conirot contro 
(Zoi) 
The complete expanded 6DOF EOM, including physical parameters for the NPS Phoenix 
AUV, are outlined in Appendix A. 


I. DEVELOPMENT OF LONGITUDINAL SURGE MODEL 


Restricting the motion of the vehicle to surge only, the significant 
motions/quantities that must be incorporated to effectively model the vehicle in the 
longitudinal direction are, the surge velocity u,, and the global position X. This 
restriction simplifies the twelve previously developed equations to a system of two 
non-linear equations of motion. Based on the NPS PHOENIX AUV equations of motion, 


see Appendix A, this reduced set which models longitudinal surge motion is, 


u(* : 20, zs prop 
§ ; (2.92) 





(m—X,)u, +X 








ulu| *r 
X=u,+U, 
This set of equations, with a slight modification, is also valid for modeling relative 


motion. By representing all variables in body fixed quantities, the set of equations for 


W-B)\,. 
a = ore 
g , (2.93) 


where the position x, is measured relative to the vehicle fixed frame. 


relative positioning is, 


(m- Ne) tee 











ulu| 


X=u,+u, 


The purpose of this longitudinal surge model is to allow for the development of a 
surge controller that will allow a vehicle to hold position in the presence of waves. To 
complete this model, the propulsion system dynamics must be included. Therefore, the 
system given in Equation 2.93 must be augmented with Equation 2.76. Augmenting 
Equation 2.93 with the propeller force equation, and assuming a neutrally buoyant 
vehicle, the set of equations to be used as the basis of the longitudinal surge controller 


becomes, 
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X=U, tu, 
u,=QUu, 


F Say, +f ut) + Ald 
T T T 








ate (2.94) 


Uu, 


where the parameters a, f, y and t must be determined through system identification. 
The process by which these parameters are determined will be presented in Chapter IV. 
Note, it should be pointed out, that in the form used in Equation 2.94, F becomes a 


generalized force with units the same as an acceleration rather than a direct thrust value. 


J. SUMMARY 


This chapter has presented the kinematic and dynamic relationships used in 
modeling a small underwater vehicle operating in a shallow water wave environment. It 
discussed the various external forces and moments that act on an AUV. The ener 
concluded with the development of the one degree of freedom (1DOF) surge model that 
will be used as basis for a controller that will allow an AUV to station-keep in the 


presence of waves. 
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IH. DISTURBANCE ANALYSIS 


A. INTRODUCTION 


In general, the disturbances that act on an underwater vehicle can be placed in 
three categories and described as follows: 
e Additive disturbances are external forces and moments which act additively 


on the vehicle. By including their effects, the total description of the vehicle 
model is extended by additional states (e.g., current, waves, and wind). 


e Multiplicative disturbances affect the dynamics of the system (e.g., the depth 
of the water, load conditions, trim, and speed changes). These disturbances 
can be regarded as time variant. 


e Measurement disturbances are due to incorrect measurements (e.g., noise on 

the vehicle's sensors). 

For this dissertation, only the additive disturbances will be taken into account. 
From the class of disturbances causing additive effects on an underwater vehicle in 
shallow water, only waves and current will be considered since they are most dominant. 
The ability to control a vehicle is known to be significantly affected by its environment. 
Since the modeling of external disturbances, especially waves, is rather complicated, 
many attempts to design control systems have suffered. 

A general assumption that is used in the modeling of AUVs is that forces and 
moments, which are added to the “‘calm sea’* model, can model environmental induced 
disturbances. This procedure was outlined in the previous chapter. This method, using 
the principle of superposition, is a good approximation for most marine control 
applications; however, it should be noted that for large general maneuvers it is not 
expected to be valid. 

This chapter will begin with a discussion of the stochastic nature of sea waves, 
including a description of several empirical relationships that can be used to represent the 
spectral content of a wave field. This will be followed by an overview of state space 
representations and recursive modeling of wave disturbances with specific application to 


control design. Finally, a methodology for the use of empirically derived spectral 
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relationships as well as measured wave elevation time series in distributed simulations 


will be outlined. 


B. STATISTICAL DESCRIPTION OF SEA WAVES 


The proper characterization of the real sea surface is difficult to obtain. With 
respect to the design of a control system, it is sufficient to assume a simplified description 
of the sea by considering only unidirectional linear waves. Based on this assumption, and 
on the superposition principle, simplified models can be determined. Concerning the 
modeling of waves, there exist two typical approaches: regular waves and irregular 
waves. 

Regular waves can be represented as a simple two-dimensional, sinusoidal wave 
train over an infinite water surface with infinite depth. This interpretation is based on the 
empirical observation, that the motions generated by waves have strong periodic 
components, see [Kallstrom 1979] for details. The characterization of the sea level as a 
train of regular waves, is however, an approximation, which is not necessarily accurate. 

Irregular waves allow the stochastic nature of waves to be taken into account. 
According to this method, the sea level can be modeled as either a superposition of a 
large number of regular waves of different amplitudes, frequencies, phase angles and 
directions of propagation, or as a narrow band stochastic process. In the case where the 
sea level is regarded as a stochastic process, the spectral and probability density 
description should be available, thus the model for the variation of the wave elevation can 
be determined from its spectral density function. The irregular waves are considered 
probabilistically with respect to amplitude and wavelength. Since the origin of waves is 
usually due mainly to the wind, the frequency and the steady state amplitude of the waves 
depend on the mean value of the wind speed. In this and the following sections, the 
stochastic characteristics of waves are considered. 

The wave elevation of a long-crested irregular sea propagating in the positive 


x-direction can be written as the sum of a large number of wave components represented 


by, 
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N 
(x,t) =>) A, sin(k;x-@,t+9;) cx 


i=] 
where, A; is the wave amplitude, @ is the wave frequency, @; is a random phase angle and 
27 


—————= OZ 
=f 3.2) 


is the wave number, with L; as the wavelength. The wave amplitude can be expressed as 
a wave spectrum S(q@), by 

A? =2S(@,)A@ (3.3) 
where, Aq@ is a constant difference between successive frequencies. The instantaneous 


wave elevation has a Gaussian distribution with a zero-mean and a variance o” defined 


as 
o =| S(@)da. (3.4) 
1. Wave Spectral Densities 


Several formulations of wave spectral densities have been proposed. The four 


spectra commonly encountered in practice are: 


e The Amplitude Spectrum. In this spectrum, the ordinates of the spectral 
density are proportional to the amplitude squared of the component waves. 
The area under the spectrum curve, Sy4, is proportional to twice the average 
energy of the record. 


e The Energy Spectrum. The ordinates of the spectral density are proportional 
to half the amplitude squared of the component waves. The area under the 
spectrum curve, Sz, is proportional to the average energy of the record and 
equals S,/2. 


e The Height Spectrum. The ordinates of the spectral density are proportional 
to the height squared of the component waves. The area under the spectrum 
curve, Sy, equals 4S,. 


e The Double Height Spectrum. The ordinates of the spectral density are 
proportional to twice the height squared of the component waves. The area 
under the spectrum curve, S2x, is therefore equal to 8S,. 


Graphical representations of these spectra are shown schematically in Figure 3.1. 
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Ze Statistical Description of Wave Amplitudes 


The probability density function of wave amplitudes with a narrow band spectrum 


can be expressed by the Rayleigh distribution, 
Dig i: 
P(r) = =e (3.5) 


where p(r)dr is the probability that a wave amplitude (r) lies between r and r+dr, and Fr’ 
is the mean square value of the wave amplitudes in the record, [Longuet-Higgins 1953]. 
It has been shown through the use of histograms that actual wave amplitudes closely 
follow this theoretical distribution, therefore, with the use of Equation 3.5, some 


quantitative statistical results may be formed. 
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Figure 3.1 Wave Spectral Density Comparisons [Berteaux 1976] 
a) Most Probable Wave Amplitude 
The most probable wave amplitude is the value of r for which 
d 
— p(r)=0. (3.6) 
dr 


Differentiating the probability density function, Equation 3.5, and equating it to zero as 
indicated in Equation 3.6, yields an expression for the most probable wave amplitude 7,,, 


aS 
r. =0.707V7° (3.7) 


where V7r~ is the root mean square value of the wave amplitudes in the record. 
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b) Mean Amplitude 


The fraction f (0 < fs 1) of wave amplitudes larger than a given amplitude 


ro 1s represented by 

f=] p(rddr. (3.8) 
The average of the f highest amplitudes can be found from the integration 

r, =< [yt ar . (3.9) 


_ The mean amplitude of all the waves in the record is then obtained when f = 1 and r, = 0 


and is given by the first moment 


r= rp(r)dr. (3.10) 


Table 3.1 shows the integration results for several values of f. 


Fraction of Largest Mean Values 
Amplitudes : (F, felr a ) 
Considered 

0.01 22350 
Gal 1.800 
0.333 1.416 
0.5 e236 
1.0 0.886 


Table 3.1 Wave Amplitude Means 


C) Maximum Expected Wave Amplitude 


The expectation of the largest amplitude in a sample of N waves is found 
from the first moment of the probability distribution of the maximum amplitudes, nx. 


Results obtained from this computation are summarized in Table 3.2. 


a5 


5. Empirical Formulation of Sea Spectra 


Several empirical formulas, based on the analysis of many wave records have 
been proposed to express the spectral density of the energy spectrum as a function of the 
wave frequency. The two most general cases use either wind speed, or significant wave 
height and significant wave period in the empirical formulas. The significant wave 
height and significant wave period, for most ocean engineering applications, is defined as 
the mean of the one-third highest waves and the mean of the wave periods associated 


with the one-third highest waves, respectively. 


Number of Maximum Wave 
Waves Amplitudes 
N (Gi! VF? ) 
50 ZZ 
100 2.28 
500 2.61 
1,000 2.78 
10,000 Sal3 
100,000 3.47 


Table 3.2 Expected Maximum Amplitudes [Longuet-Higgins 1953] 


In the first case, the spectral density S(q@) is of the form 
A 4,44 
S(@)=— eC" en (3.5) 
@ 


where A and B are empirical constants, and V is the speed of the wind. 
In the second case S$(@) is given by 


INGLE ie 
S(@) = - xe Ee 


(3.6) 


where A and B are again empirical constants, 7; is the significant period and H, the 


significant wave height. Table 3.3 presents a useful compilation of wave heights and 
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wave periods as a function of sea states and wind speed that may be used in the following 


formulas. 


nye) Wind Wave Height Wave Period 
State Velocity (feet) (sec) 
Cen) 


Average — Significant TAT I lice 


eh Batak aoe or 
0.10 Si |r, 
= 





1.4 





ea A 


10 13 3 AD il iS Te 
24 7.9 i 16 Beale 6.8 
Ans a2 13 17 3.8-13.6 | 9.9 7.0 
26 9.6 15 20 4.0-14.5 | 10.5 7.4 
28 11 18 23 A515 Sele 7.9 


Table 3.3 Sea State and Wave Parameter Comparison [Berteaux 1976] 


a) Pierson-Moskowitz (PM) Formula 


A frequently used one parameter description of S(@) for a fully developed 


sea, resulting from extended wind with unlimited fetch, is the PM spectrum. Using 
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Equation 3.5 as a basis, the single parameter used to describe the spectrum is the wind 
velocity V, with the empirical constants A and B having the values 0.0081g” and 0.749%, 
respectively. The units associated with the gravitational constant g, and the wind velocity 


V must be the same for dimensional consistency. 
b) International Towing Tank Conference (ITTC) Formula 


Using statistical properties the relationship between the one-third 


significant wave height, 1/3, and the wind speed, V, for the PM spectrum is 
V 2 
igi — (I) 0s —— (3.9 
§ 
therefore, the formula for the PM spectrum with the significant wave height as the single 


parameter is 
oe yur 
S(@)=—<e Cea) (3.8) 


In this form, the empirical constants A and B have the values 0.0081g * and 0.0324 ”, 
and the resulting one-parameter expression that uses significant wave height 1s referred to 
as the ITTC formula. 

Since this one-parameter formula describes a fully developed sea resulting 
from conditions that are rarely encountered (extended wind with unlimited fetch), the PM 
spectrum should be viewed as an asymptotic form. To overcome the limitations 
associated with the one-parameter spectral family, a two-parameter family, given by 


Equation 3.6, can be used. 
Cc) Bretschneider Formula 


The Bretschneider spectrum is a general form, and represents experimental 
data very well. The formula for the Bretschneider spectrum is 


2 
4200H 5 050K T} a4) 
Te aye : 


s 


S(Q@) = (3.9) 


where H, and T; are the significant wave height and significant wave period, respectively. 
This expression can be used to represent developing, fully developed and decaying seas 


depending on the value of 7; chosen, [Lewis 1989]. 
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d) International Ship Structure Congress (ISSC) Formula 


The ISSC spectrum is similar in form to the Bretschneider spectrum, with 
the exception that it is based on a mean frequency corresponding to the spectrum's center 
of area. This spectrum is represented by 

Ip eae errr 
S(@) == =e ela (3.10) 
where H, and T are the significant wave height and mean wave period, respectively. A 
comparison of the ITTC, ISSC and Bretschneider spectrum for a sea state 3 condition is 
shown in Figure 3.2. 


1.8 


ITTC spectrum 
16 = .. ISSC spectrum 
oe Bretschneider spectrum 


0.4 


0.2 





@ (rad/s) 


Figure 3.2 Comparison of Empirical Spectra 


e) Other Spectral Representations 


There are many other empirical relationships that have been used to model 
spectral characteristics of the sea, where each model is useful for a given area or sea 


condition that has specific characteristics. Examples of these are the Joint North Sea 


a) 


Wave Project (JONSWAP) spectrum [Hasselmann 1973], the Ochi Six-Parameter Wave 
spectrum [Ochi 1076], the Wallops spectrum [Huang 1981] and the Generalized spectrum 
[Liu 1983]. The Wallops spectrum and the Generalized spectrum with variable 
exponents are adaptable for both deep and shallow water applications. 

It is common to use the recommended sea spectra from the ITTC and ISSC in ocean 
engineering design. For open sea conditions, the PM spectrum is recommended, and for 
fetch limited conditions either the Bretschneider, Ochi or JONSWAP spectrum is 


available. 


c. LINEAR REPRESENTATION OF SEA WAVES 
1. Spectral Approximations 


To properly estimate and cancel the wave induced disturbances acting on an 
underwater vehicle, some type of disturbance model is necessary. The necessity arises 
from the need to either embed the disturbance model in the vehicle's controller, or to use 
the disturbance model in a State estimator. It is known that a linear, Gauss-Markov 
stochastic process may be generated by sending white noise through an appropriate 
transfer function, where white noise is to mean a random process whose power spectral 
density is constant over the whole spectrum. Therefore, a linear approximation to the 
spectral representations in the previous sections can be obtained by sending a random 
white-noise signal through a second order filter, [Spanos 1981]. This process can be 
written as 

y(s) = h(s)q(s) Gi 


with the linear approximation to the desired spectrum represented as 


®, (a) =|h(ja)|®,, =|Ajo)|’. (3.12) 


In Equation 3.11, y(s) is the wave amplitude, g(s) is a white noise source with power 


spectrum 


®,, =1.0 (3.13) 


and h(s) is a second order transfer function of the form 
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2k(s/@,) 


~ 1422(s/@,)+(s/@,)° ee 


h(s) 


Using equations 3.11-3.14 and substituting s = j@ the approximate power 
spectrum can be given by 


4k*(@/o,)’ 


ee ee (3 
(l-(@/@,)*)’ +46°(@/0,)’ oi) 


® ,, (@) = |h(jo)|’ ® ,, (@) = 


The power spectrum of this filter, ®,,(@, has zero energy for zero frequency and the 
maximum value of ®,,(@ occurs at @ = @, for small values of ¢ This is desirable 
considering the shapes of the spectra displayed in Figure 3.2. The parameters ¢ and k in 
Equation 3.15, are found by minimizing the performance index 
J=| > (®,, (ja)-S(@)y do. (3.16) 
ja State Space Formulation 


A linear state space model can be derived from equations 3.11 and 3.14. By 


defining the states as 


=[ y(r)dr 
ai (3.17) 
xX,=y 
the system states can be represented by the set of equations 
X, 0 1 x v 
a 2 a q 
Ke -@, ~-2¢6@, |x, | | 2ka, 
: (3.18) 


ro le 


The transfer function in Equation 3.14, is useful for representing the sea surface 
elevation, but it cannot be used to generate the wave velocity. This is seen by taking the 
limit 

lim sy(s), (3.19) 
as s tends to infinity. The result of this calculation is a constant, 2kq@,, that cannot 
represent the oscillatory wave velocity. [Saelid 1983] has shown that using the proper 


transfer function 


ao 


y_ 2¢0(s/@,)s (3.20) 

q (1+2¢(s/@,)+(s/a,)*) . 
can solve this problem and [Riedel 1997} has demonstrated through the use of linear 
prediction theory that an eighth-order transfer function of the form 
Ui, -| 260 0)s 


q 1+2¢(s/@,)+(s/,)° Ce 


allows extremely accurate matching of the target spectrum enabling fluid velocity 
prediction as much as one period ahead. 


Based on Equation 3.20, a generalized transfer function may be written in the 


form 
—— (3.22) 
q (atbs+s*)’ 
which implies 
(s* +2bs° +(2a+b*)s* +2abs+a*)u ,(s) =kq(s). (3.23) 


[Astrom 1989] has shown that the parameters a, b and k, representing a PM spectrum 
based on sea state 3 conditions, can be found using Equation 3.16 to be a = k = 1 and 
b=2. In this form, Equation 3.23 becomes an all-pole filter thus avoiding the problem of 
base period repetition of wave records, [Riedel 1997]. 


25 Spectral Modifications for a Moving Vehicle 


The spectrum that the vehicle “‘sees” while moving through a wave field with 
some forward velocity is not the same as that which a still vehicle would encounter. The 
actual spectrum that the vehicle encounters is a function of the vehicle's forward speed U 
and its heading angle relative to the propagation direction of the sea waves, f. The 
definition of the encounter angle / is given by 

B=n-(y-y), (3.24) 
where y is the direction from which the waves are propagating referenced to the inertial 
reference frame, and wis the heading angle of the vehicle [Lewis 1988b], see Figure 3.3. 


The frequency modification that must occur to the disturbance spectrum is represented by 
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the Doppler equation, (Equation 2.56), given in the previous chapter. The Doppler 


equation, in an alternate more useful form is 


Z 


o, ~@-—-U cos B. (3.25) 
§ 

Depending on the encounter angle, specific terms have been given to the 
orientation of the sea with reference to the vehicle, this is shown in Figure 3.4. As an 
example, consider the case where the encounter angle is zero, namely the vehicle is 
moving in the direction of wave propagation, this is referred to as a ‘‘following sea’*. In 
this case, as equation 3.25 indicates, the frequency of encounter, @, is less than the actual 


wave frequency. For this case, it is interesting to note that the encounter frequency can 


be negative for large values of forward velocity U. 


wave 
direction 





B 





Figure 3.3 Incident Wave Directions 


The methods for spectral approximation and state space realization of sea waves 
presented in this section are useful in simulation studies when the target spectrum is 


known, including the frequency shift. However, if the target spectrum is not known, 
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which is typically the case with a deployed vehicle, alternate techniques must be used. 
The following section presents one such technique which will allow a vehicle to estimate 


the encounter spectrum on-line. 


B=902 
Beam sea 







p= 45° 
Quartering sea 


= [350 
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Figure 3.4 Encounter Angle 


D. WAVE MODELING USING RECURSIVE METHODS 


A typical experiment in system identification consists in recording a set of 
input/output data and fitting a parametric model. In discrete time, the attempt is to fit a 
Linear Difference Equation (LDE) model of the form 

yO) +a yi—1)+...4a,yG —-n) = bree + OU an) eh), (3.26) 
with t:[0,°0] denoting the integer discrete time index, and at) an error term which 
accounts for the fact that the data never matches the model exactly. The problem 1s to 
estimate the parameter vector 

Ge Nomaa DaseeDal (3.27) 
from a set of data. This section will present the general procedure to estimate the 
parameters associated with a discrete linear model. 

Equation 3.26 may be written in regression form as 


y(t) = O(t-1)'0 + e(t) (3.28) 
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with 
o(t-—1) = [— y¢—-D),....-y(t—n),u(t—-D),...,.u(t-—n) (3.29) 
being a sliding window of input and output data. Now the problem is to determine a 
technique to compute an estimate of the parameter vector @ on the basis of the data set 
ZZ = {u(0),. .. U(N), y(Q),..., \ON }, if it is assumed that N data points are collected. 
If the problem is cast in a probabilistic framework, a probability density for the 
data set Z given the model @ can be written as 
P(Z@) = Pr(€(0),...,€(N)) (3.30) 
with €(t)= y(t)-¢’(t-1)@. If we assume the disturbance term sequence to be 


Gaussian and white, then the density on the nght hand side becomes 
N 
Pr(€(0),...,E(N)) = | | Pr(e(a)) (3.31) 
t=0 


with 
e? | 


oom (3.32) 





ae a 


As a consequence it can be seen that 


5 DI 9-0" D9 
Pr(Z/6) = Ce | (3.33) 


Minimizing the summation in the exponential can maximize this probability. So 


if the estimate of the parameters is defined as the vector that minimizes the probability, as 
6 = arg min, Pr(Z{6) (3.34) 


then it can be computed as a least squares solution 
uM 2 
6 =argmin, > |y(t)- 97 (2)6 (3.35) 
t=] 


The solution can be obtained using standard techniques, by writing 6 as the least 


squares solution of the system of equations 
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y() o(1)" 


a = g2y" 0 (3.36) 
y(M)} | ecm)’ 
or, compactly 
y=OD0. (3:39) 
The solution given by the pseudoinverse becomes 
6=(0'®)'o'y (3.38) 
which can be written as 
ee ie 
i dome] ZF Eeor0 (3.39) 


This is true provided the error sequence &t) is white and Gaussian. If the error is not 
white, then the estimate of the parameter vector @is biased. 
Now that we have set the foundation, the problem that arises is how do we 


estimate the parameters associated with a transfer function that will properly represent a 


model of the seaway in a recursive fashion. If O(t)is the estimate of the parameters at 


time t, the goal is to compute the successive estimate 6(t + 1) by updating 6(t) using the 
latest observations. It turns out that this problem can be put in a very nice framework that 
makes use of the considerations on the Kalman Filter. More on this approach will be 
discussed in Chapter IV. 
The auto-regressive (AR) model, with numerator equal to one, can be written as 
y(t) = @(t-1)' 6+ e(t) (3.40) 
with e(t) a white noise sequence. If it is assumed that @ is constant, Equation 3.40 can be 
written in state space form, where the state is the parameter vector itself, as 
O(t+1) = A(t) 


. | (3.41) 
y(t) = Gt — 1)’ A(t) + e(t) 


This form is just a particular case of the stochastic state space model with A being the 


identity matrix, B=0 and C=¢(t)’. Since this is a stochastic process, the Kalman Filter 


approach for the estimation of &t) is used, which leads to 
P(t)g(t -1) 

A +0(t-1)’ P(HeOt-1 

P()o(t-Nge@-1)" Pi) 

M+ G(t-1)" P(t)o(t-1) 


(y(t) - 8(2)" ot -1)) 
, (3.42) 


6(t+1) =O(t)+ 
P(t+)l=P(t)- 


with A” = Elec)’ f Clearly, the parameter A is never known, and the need to know it 


can be eliminated by proper normalization. Dividing the numerator and denominator of 
Equation 3.42 by A’ the following algorithm is obtained; 
P(t)o(t —1) 
1+@(t-1)" P(t)e(t-1) 
P(t)o(t- Gt)" P(t) 
1+@(t-1)’ P(OCt-1) 


which can be easily verified by setting P(t)=P(t)/2*. This algorithm is called the 


(y(t)- A(1) @(t-D) 
: (3.43) 


O(t+1)=6(t)+ 


P(t+1)=P(t)- 


Recursive Least Squares (RLS) algorithm, [Ljung 1987]. 

Using the above developed RLS algorithm the parameters of an AR model of the 
form 

w(t) +d,w(t—1)+...+d,w(t—N) = e(t), (3.44) 

with w(t) representing the sea surface elevation due to wave action, and e(t), a zero mean, 
white noise otherwise know as the innovation, may be developed. The main advantage of 
the AR model is the fact that the parameter estimation is a linear operation. In fact, if the 
numerator pertinent to the noise term is not one, then the error state model, Equation 
3.41, is not white and we cannot apply Kalman Filtering techniques directly. In the 
/ general case, the problem is nonlinear in the parameters and the Extended Kalman Filter 
must be applied. 

The problem that now arises is to determine the order (Y)) of the model to properly 
reflect the actual frequency spectrum of the time series w(t) while keeping the complexity 


of the model at a minimum. By computing the covariance of the innovation as the order 
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of the model is increased it can be shown that there is a value to which the covariance 
converges. Using surface wave elevation data obtained in Monterey Bay, from a 
Waverider measurement buoy, an Auto Regressive model of various orders was 
determined and the covariance of the innovation compared. The results of this analysis 
are shown is Figure 3.5. 

As can be seen, the covariance begins to flatten out around an eighth order model. 
This would have one believe that the correct order model to choose would be an eighth 
order model. However, as shown in Figure 3.6, an eighth order model fails to accurately 
reflect the actual spectrum. As the order of the model is increased, the error in the 
matching of the actual spectrum decreases, however, not until a 100" order AR model is 
computed, is the desired spectrum actually realized, see Figure 3.7. The issues associated 
with using a model order this high include parameter corruption due to noise causing 


inaccurate identification. 
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Figure 3.5 Comparison Of Innovation Covariance To AR Model Order 


There are those that will argue that the energy associated with the low and high 
frequency modes is minimal, and that proper modeling of those modes is not important. 


However, as was discussed in Chapter II, it is the low frequency modes of the wave train 
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that are considered shallow water, non-dispersive waves. It is these modes that will have 


the most effect on the submerged vehicle trying to maintain station. 
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Figure 3.6 8'" Order AR Model — Fit to Monterey Bay Data 
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Figure 3.7 100" Order AR Model- Fit to Monterey Bay Data 
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Therefore, as in any design, trade-offs must be made based on analysis, as to how 
accurate the wave disturbance model must be compared to the resulting controller 
complexity. 

Using a model order much greater than that of an eighth order approximation 
creates difficulties in real-time embedded processes currently used on small underwater 
vehicles. However, using a stable, reduced order disturbance model in an estimator 
where measurements are available to correct the model errors provides a very good 
means of tracking and compensating for the disturbance, [Riedel 1998a]. Further 
information regarding the estimation and accuracy of these models will be discussed in 


Chapter V. 


E. APPLICATION TO DISTRIBUTED SIMULATIONS 


This section will present a method by which realistic wave data as well as 
empirical relationships may be implemented into a distributed simulation. The purpose 
of this approach is to further the development of simulation capabilities for control 
system design, multiple vehicle coordination as well as mission planning. 

The simplest method of obtaining information about wave disturbances in a 
particular operating area is from a wave buoy. However, to use this information to 
simulate the disturbance forces and moments acting on a submerged vehicle, it 1s 
necessary to transform the wave elevation record into a water particle velocity record at 
the vehicle operating depth. There are two approaches to this problem. The first uses 
spectral analysis while the second uses Fourier analysis. 

Using the spectral analysis approach, the procedure is to first compute the power 


spectral density, Py, @), of the surface elevation, then modify the PSD by the appropriate 
depth related transfer function, \H (OeZ 1, and finally convert the new spectral density 


back to the time domain for replications of the subsurface water particle velocities. 
Although the resulting time series reflects the proper magnitude of the water particle 
velocities, the disadvantage to this method is that the phase information 1s lost, and the 


resulting subsurface velocity time series do not accurately reflect the motion caused by 
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the surface elevation times series. It is critical to ensure that the phase relationship 
between the horizontal and vertical wave induced velocities and accelerations are correct 
in order to provide realistic vehicle motion in simulation since, as Chapter JI outlined, the 
6DOF EOM are coupled. 

Through the use of the Fourier analysis method, the shortcomings of the spectral 
analysis are overcome and the disturbance phase relationship is maintained. In this 
method, a FFT of the surface wave record is taken, then the Fourier coefficients are 
multiplied by the appropriate value of the modifier for each frequency component. Once 
this has been completed, an inverse FFT is performed. The resulting disturbance records 
now reflect both the proper amplitude and phase relations. Figures 3.8-3.11 depict this 
Fourier analysis translation procedure for a wave elevation time series recorded in 


Monterey Bay, CA. 
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Figure 3.8 Wave Elevation Time Series, Monterey Bay April 1998 
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Figure 3.9 Surface Elevation PSD, ®,,(/) 
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Figure 3.11 Horizontal Water Particle Velocity PSD, ®y,(f), (H=45 m, Z=22.5 m) 


The elevation data, displayed in Figure 3.8, was recorded by a Waverider® buoy 
in 45 meters (150 feet) of water. The resulting velocity record, Figure 3.10, is from a 
transformation of the elevation record for an operating depth of 22.8 meters (75 feet). 
Referring to the PSD of the transformed velocity record, Figure 3.11, it is interesting to 
note that in addition to the dominate peak frequency at 0.12 Hz , there are two additional 
lower frequency components at 0.05 Hz and 0.07 Hz. Recalling comments made with 
respect to shallow water waves in Chapter II, for this water depth only the frequency 
component less than 0.05 Hz may be treated as shallow water waves. 

To accurately reflect, in simulation, the wave field approximated by any of the 
empirical relationships discussed in Section B, Equation 3.3 as well as the relationships 
given in Table 2.2 are used. The stochastic nature of the sea waves is introduced in the 
simulation by including a random phase angle in the oscillatory @ term given in Table 
2.2. The proper velocity amplitudes for this wave field are determined by using the 


vehicle’s depth in the exponential modifier. The correct phasing, with respect to the 
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vehicle, is obtained by using the vehicle's global position, X, in the @ term associated 
with the water particle velocities. As a note, this approach models first order wave 
processes, however, to adequately represent the nonlinear second order wave effects, the 


Aqgin Equation 3.3 is often varied to provide frequency bands of equal energy. 


1 SUMMARY 


This chapter has presented a detailed discussion of wave disturbances. It has 
described the statistical nature of the sea and presented several empirical relationships 
that may be used for approximating the spectral properties of the ocean. Methodologies 
for identifying and employing state space models of the sea, in control applications, have 
been highlighted. Approaches for using wave information, from empirical relationships 


or real data, in distributed simulations has been discussed. 


IV. PARAMETER ESTIMATION 


A. INTRODUCTION 


As shown in Chapter II, accurate modeling of underwater vehicles has lead to the 
development of complicated equations of motion. Apart from the non-linearities inherent 
with underwater vehicle motion, the forces and moments acting on the submerged body 
are typically determined by a combination of theoretical and experimental results. 
[Healey 1992 and 1993] has suggested the use of three, independent, decoupled, 
equations of motion sets, which model speed, steering and diving control. The use of 
these simpler models, which describe only particular vehicle dynamics, is useful in 
control law design. 

Recent interest in underwater vehicle maneuvering and control in shallow water 
has generated the need for a greater understanding of vehicle dynamics in this regime. 
Specifically, improved vehicle models foster the development of sophisticated control 
architectures, which produce the high degree of autonomy necessary to allow vehicles to 
maintain acceptable performance in an ocean environment. Critical to the solution, since 
dynamic positioning of an underwater implies a nonlinear response, is the use of an 
adequate input-output mapping of the vehicle dynamics. 

In this chapter, a method for identifying the decoupled longitudinal surge motion 
dynamic parameters is presented. The identification is based on post-process Kalman 
filtering of data obtained from in-water vehicle experiments. Identification of the 
parameters associated with the surge equations of motion is performed, and a comparison 
between experimental data from in water measurements with the Phoenix AUV, and 
simulated results is conducted. Lastly, the nonlinear function that relates the commanded 
propeller speed to the required propulsion motor voltage is determined. This function is 
critical in the implementation of a real-time surge controller that will allow a vehicle to 


maintain position while disturbed by waves. 
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B. ESTIMATION METHODOLOGY 


Parameter estimation has been called a “can of worms” by Astrom 
[Astrom 1983], in which he refers to the difficulty of making theoretically sound 
methodology work with real data. With this said, many different techniques have been 
employed in the area of system identification or parameter estimation, see [Lung 1987, 
Gelb 1968 and Astrom 1989] for examples, to attack this daunting task. In this work, a 
recursive Kalman filter approach was chosen since this technique is suitable for real-time 
implementation. The Kalman filter method is similar to weighted least squares 
algorithms [Ljung 1987], and allows for the incorporation of system and measurement 
errors. 

In the Kalman filter algorithm, it is assumed that the parameter model is based on 
a nominally constant parameters set, where the state vector is the parameter vector 0, and 
the dynamics are described in discrete time by 

6 (k+1)=0 (k)+Iw | 4.1) 
2(k +1) =h(k)O (k)+v 
The system noise, w, and the measurement noise, v, are considered to be zero-mean, 
white noise sequences with associated covariance matrices Q and R, respectively. The 
recursive Kalman filter estimation equations as given by [Gelb 1974] are 
6. Demo 
ee ere een 


kik-t 


K, = Py Ee, es +RJ]~ ’ (4.2) 


kik Oem ta [Z, —h, 8, ;,-4] 
P,,, =U-K,A, VP 


+ Gy oOle 


kik=-1 


5) 


where K is a time varying optimal gain that produces a least squares solution for the 
parameter estimate, and P is the parameter estimation error covariance. These 
expressions are equivalent to the formulation by [Ljung 1987], where the ‘‘forgetting 


factor’ ‘ A, is related to the noise covariance R. This recursive algorithm is expressed as 
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6(t) = A(t —1)+Litje(t) 
e(t) = (z(t) h(t -)] 
L(t) = P(t)®(t) = P(t -NOMADI+@" (NPC -NOew@]'> 4-3) 


a 
P(t) = 1 P(t-1)- P(t = (t)P(t—1) 
A(t) A(t)E + O° (t)P(t- YO) 
where the standard Recursive Least Squares (RLS) algorithm is obtained for the special 


case of A= 1. 


C, PARAMETER IDENTIFICATION MODEL 


The surge model developed in Chapter II is in continuous time, however, 
real-time control and estimation is done in discrete time, therefore the set of equations 
given by Equation 2.91 must be converted to a digital form. To produce a digital set of 
equations, a standard Euler discretization can be used. Assuming a sampling period 7, 
the discrete system of equations, disregarding the kinematic relation, for the surge model 
becomes 


u,(k +1) = (au, (k) 





u, (k)| + F(k))T' +u,(k) 
(4.4) 


F(k+l= Srittu, (k)|In(k)| +E n(kace pr” + F(k) 

The model presented in Equation 4.4 contains four unknown parameters, a, f, vy 

and 7, that must be determined by experimental means. To properly identify the dynamic 
parameters associated with this mode, the model must be expressed in terms of the 
measurement variables only. The measured input-output data channels that are available 
for this identification are the relative longitudinal velocity, u,(t), and the propeller 
revolutions, n(t). Since the propeller thrust value F, cannot be measured, the two 
first-order equations must be combined into a second-order equation containing only 


these two measured variables. Defining the change of variables 
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poe 
c 

aU (4.5) 
T 

pen 
- 


and combining the two first-order equations in Equation 4.4, a second-order discrete 
model in u, can be formed, 
u,(kK+2)=2u_(k+1)—u,(k) 
+Tatlu,(k+1u,(k+1|—-u, (lu, GO] 
+Tdlu, (k)—u,(k+D]+T?balu, (ku, Gol] 
+Tclu, (K)\nGo|]+T7[ndo|nce)]] 








(4.6) 





Unfortunately, the model in Equation 4.6 is nonlinear in parameters making it 
difficult to apply standard system identification techniques. However, the nonlinear 
model can be transformed to be linear in its parameters and variables, thereby allowing 


the use of well established estimation tools, by defining an additional change of variables 











as, 
C.=Te h, =u, (k+)lu, (k+)|-u, (ku, (| 
C, =Tb h, =u,(k)—u,(k +1) 
Cre ie h, =u, (k)|u, (k)| (4.7) 
Cae h, =u, (k)\n(k) 
C2720 h; =n(k)|n(k)| 
and 


Z=uUu,(k+2)—2u_(k +1) +u,(k) 
h=|, bh hy hy hs) (4.8) 
b=|C, 16m GC Vemelt 
With these definitions, the model can now be written in matrix notation as 
z=hé. (4.9) 
The parameter vector, 6, contains five coefficients where only four parameters are of 


interest. The additional coefficient, C3, is a result of the nonlinearity associated with the 
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original model represented by Equation 4.9. The parameter vector, @, will be estimated 
in a least squares sense using Equation 4.9, and in this case, one extra degree of freedom 
is present which may cause a slight decrease in the accuracy of the other parameters. 

In theory, the additional coefficient defines an implicit relationship between the 
other parameters of the model; in reality however, this additional coefficient lumps 
together any modeling errors and therefore it is not constrained thereby allowing 
convergence of the parameters in a least squares sense. The model parameters @, f, yand 


T, are therefore determined from the four coefficients C;, Cs, Cz and C2, respectively. 


D. SYSTEM IDENTIFCATION 
1. Input Signal Design 


Prior to estimating the parameters associated with the surge model, data must be 
obtained for use in the estimation filter. This data must contain measured control input 
and response output variables from a “‘persistent’* excitation. It is a well-documented 
theorem that to ensure a unique, unbiased, least squares estimate, the system or plant 
must be persistently excited, [Astrom 1989]. In addition, the system identification should 
take place using open loop control, if possible, so that controller dynamics do not effect 
the results, [Ljung 1987]. 

Since the purpose of this piece of work is to determine the parameters of the surge 
dynamics model that will be the basis for a surge controller, the input data must 
persistently excite the vehicle over the expected frequency range of the surge velocity 
disturbance. For shallow water applications, the period of the surge disturbance that an 
underwater vehicle may encounter can range from approximately 4 to 40 seconds. It is 
necessary therefore, that the propeller revolution input to the vehicle, for identification 
purposes, also contains this frequency content. By selecting a square wave of various 
periods the control input was designed that contained the desired frequency components. 
A portion of the time series used for the control input is shown in Figure 4.1, and the 


frequency content of this input signal is depicted in Figure 4.2. 
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Figure 4.1 Control Input Time Series 
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Figure 4.2 Control Input Frequency Content 
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As seen in Figure 4.2, the frequency components of this input signal range from 
0.02 Hz to 0.35 Hz. The dominant frequency was designed to be around 0.05 Hz to 
0.08 Hz,or a 12.5 to 20 second period. The purpose of designing the input signal in this 
manner is to ensure that the parameters were weighted in the range of the dominant surge 
period. With the input signal properly designed the system identification data runs can 
proceed. Note that considering the input to be n(t) or n(t)|n(t)| makes no difference in the 
spectral content of the input signal. 


2. Data Collection 


A series of four in-water experiments were conducted in the Monterey Harbor 
Basin on March 5, 1999. The Phoenix vehicle was placed under active control in both 
heading and depth through the use of control surfaces, and the propeller RPMs were 
commanded through the use of an input file representing the input signal shown in Figure 
4.1. During these runs Phoenix carried its standard sensor suite, which includes a 
SonTek® ADV and a RDI® Navigator DVL. These sensors were used to measure the 
vehicle's response to the control input. Information concerning each of these sensors may 
be found in Appendix B. A sample of the measured input and output data obtained 
during run three is depicted in Figure 4.3. 

The input-output response shown in Figure 4.3 indicates that the vehicle is less 
efficient when operating astern. This is evident by observing that the magnitude of the 
astern velocity is significantly less than that of the forward velocity, for the same 
propeller revolutions. This decrease in efficiency is an issue that must be addressed 
during the estimation process. 

An additional item observed during data analysis is that the voltage required by 
the propulsion motors to obtain the same revolutions is different. The reason for the 
difference is attributed to the starboard propulsion train having to overcome a greater 
amount of friction. This friction is caused by misalignment that resulted from structural 
damage that was incurred by Phoenix early in its operating life. The voltage difference is 
shown in Figure 4.4. To account for the differences in propeller revolutions, an average 


input was used during the system identification. The speed from the two shafts was 
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Figure 4.4 Digital Voltage to Propulsion Motors 


averaged, and the vehicle was treated as having a single input. The parameter / is then 
determined for the assumed single shaft system, with the recognition to the fact that the 
actual parameter for each shaft on the vehicle will be //2. 

Since the output of the digital controller is a commanded voltage, this voltage 
mismatch can be accounted for in the controller software. To ensure that the propellers 
produce approximately equal thrust, the function that maps the commanded revolutions 
from the controller to delivered voltage to the motors must be determined. Once this 
function is found, it can be coded into the vehicle control computer thereby ensuring the 
propeller mismatch is minimized. Another way of ensuring that the shaft speeds are the 
same is to “‘close the loop’* around propeller speed. However, this solution adds 
additional system dynamics making the controller formulation more difficult. The 
methods used to determine this mapping will be discussed later in this chapter in 


Section E. 


S: Parameter Identification Results 


In the application of the parameter estimator, the measurement noise covariance, 
v, was set to 0.01, a constant scalar. The values of the diagonal elements, g;;, of the 
model noise covariance matrix, Q, were chosen to match the bandwidth of the input 
signal. With this choice, a tradeoff between convergence, stability and precision of the 
estimates 1s made. The following results are presented for one of the many data runs and 
selected values of gj that produced the “best” results in parameter estimation. The term 
“best” 1s a Subjective measure based on a balance between the whiteness of the residuals, 
the values of the parameters compared to parameters previously identified and a 
comparison between measured data and simulated results. The final results of the 
selected parameters from the estimation process are shown in Table 4.1, including 
Statistics on the estimation error and the residual e. 

The evolution of the parameter estimates and the diagonal elements of the 
covariance matrix, P, associated with these parameters is shown in Figures 4.5 and 4.6, 
respectively. Referring to Figure 4.5, it can be seen that the parameters do converge, but 


exhibit some slight fluctuations. The noise observed in the evolution of the covariance 
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Table 4.1 Parameter Estimation Values and Residual Statistics 
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Figure 4.5 Parameter Evolution 
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Figure 4.6 Covariance Evolution 


matrix, Figure 4.6, is caused by the rather large values of the diagonal elements of Q, 
gi = 5.0. The large values needed for these were necessitated by the additional degree of 
freedom in the filter, the uncertainty associated with the thrust reduction term and the 
large bandwidth of the input signal. 

The performance of the filter can be determined by computing the autocorrelation 
of the residuals. The autocorrelation provides a measure of how the value of a random 
variable at a time ¢, will influence its value at some future time, t+7. If the filter is 
properly tuned, the residuals should be white. For the residuals to be classified as white, 
there will be ideally zero correlation at any time shift 7 and the time series will be highly 
correlated at T= 0. Figure 4.7 indicates that the residuals are not white, but exhibit some 
correlation. The non-whiteness of the residuals indicates the lack of modeling capability. 
Since Kalman filter theory assumes that any measurement noise has a Gaussian 
distribution, the measurement noise model is corrupted by unmodeled, colored noise. 
This corruption of the measurement model caused the residuals also to display colored 


noise properties. Also, since the reduction in propeller efficiency when the vehicle 
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Operates astern is not known, additional modeling errors are introduced which effects 
filter performance. 

Despite the issues associated with the non-whiteness of the filter residuals, 
apparent satisfactory parameter identification was obtained. The parameter estimation 
was deemed satisfactory by comparing the measured vehicle response to simulated 
response results. Using the identified parameters, listed in Table 4.1, and the measured 
propeller input shown, in Figure 4.3, a simulation was conducted with the results 
presented in Figure 4.8. 

The response traces shown in Figure 4.8 indicate that the identified parameter set 
provides a reasonable predictive response. Recalling that the parameters are estimated in 
a least squares sense, the errors between the two traces are acceptable. The two 
responses are not lagged indicating good agreement in the propulsion system time 
constant and vehicle drag coefficient. The error in response amplitude is attributed to the 


uncertainty associated with the identification of the B and yparameters. 







Figure 4.7 Autocorrelation of Estimation Filter Residuals 
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Figure 4.8 Comparison of Measured and Simulated Vehicle Response. 


The uncertainty in each parameter may be estimated by using the magnitude of 
the covariance matrix diagonal corresponding to the estimated parameter. The calculated 
uncertainty for each of the parameters, based on the covariance levels, is shown in 
Table 4.1. Although the level of uncertainty with these parameters may seem excessive, 
it is not too different from the “standard” level of uncertainty associated with underwater 
vehicle parameter estimates, which typically is around +40 %. It is this variability in 
parameters, as well as other items, that produces the need for robust control law design, 


which will be discussed in the following chapter. 


E. PROPULSION SYSTEM BALANCING 


As outlined earlier in this chapter, the required propulsion motor voltage to 
produce a desired propeller revolution is different for each shaft. Under position control, 
this difference if not accounted for could cause unsatisfactory results since for position 


control the loop is closed around position error and not propulsion shaft speed. 
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In the control system architecture of the Phoenix AUV, the shaft speed is 
controlled by a motor voltage command which originates from the execution level 
computer. The particular control law that has been implemented in the vehicle calculates 
this voltage; therefore, the relationship that maps desired shaft revolutions to voltage 
must be determined. 

By fitting a curve to an overlay of data obtained during the four, system 
identification runs, the relationship between propeller “revs” and motor voltage for each 
shaft can be determined. Using the MATLAB® polyfit algorithm it was determined that 
a third order polynomial adequately mapped the two input-output relationships. The 


nonlinear functions for these relationships are, 


V iene = 0.022607n3,,, + 0.000194n2,, + 0.832982n,,,, aa 
Vieg = 0.022356n3,,, - 0.009507n2,, + 0.345086n,,, 


with the graphical results of the “goodness of fit’ shown in Figure 4.9. 


Motor Voltage (v) 


Motor Voltage (v) 





0 
Shaft Revs (rps) 


Figure 4.9 Nonlinear Shaft Revolution to Motor Voltage Function 


As a note, the polynomial curves in Figure 4.9, indicate that the starboard shaft 


will produce approximately 8.5 rps for an applied voltage of 24 volts, while the port shaft 
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produces almost 10 rps for the same applied voltage. Although there is no data displayed 
in Figure 4.9 to validate this curve at high speed, the shaft revolution values estimated by 


these curves have been observed during full power trials. 


| SUMMARY 


This chapter has presented the work performed to identify the parameters 
associated with the longitudinal dynamics of the NPS Phoenix AUV from open water 
experiments. The results show the necessity of a non-linear dynamics model and of the 
need to include a force lag and a thrust reduction term in the propulsion model. The 
results obtained are satisfactory although the filter residuals appear to be non-white due 
to limited modeling capabilities. Regardless of the level of whiteness, it must be stressed 
that the simple model identified is very useful for control law development. 

One interesting point observed during the estimation process was that taking QO as 
a diagonal matrix with equal values produced better results in the sense that the residuals 
‘were minimized and the filter produces stable parameter evolution. This reflected a 
priori, equal certainty in each initial parameter estimate. Moreover, the interrelationships 
between the parameters implied a uniform characterization of the modeling noise in the 
filter. The tuning of the estimation filter proved to be an important step towards 
obtaining a good model. 

Finally, the methodology presented here can be extended to other types of vehicle 
motions and represents a basis for the development of models and identification 


techniques that can be applied to other coupled or decoupled motion models. 
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V. DISTURBANCE REJECTION THEORY 


A. INTRODUCTION 


Many different techniques have been applied to the disturbance rejection problem 
in the past. For years, passive techniques employing improved electro-mechanical design 
were used. With current advancements in the computing industry, adaptive methods have 
become a popular means of active control. Recent research in the area of nonlinear 
control has shown that Variable Structure or Sliding Mode Control provides a very robust 
method of disturbance compensation. 

This chapter will begin with an overview of disturbance rejection theory 
developed to date. This will provide the reader a foundation from which to see the 
extensions made here in regards to the development of disturbance rejection techniques 
for small underwater vehicles. 

| Next, a series of three case studies will be presented. These case studies will 
demonstrate the major approaches available in the design of disturbance rejection 
compensators, including the effect that the spectral content of the input disturbance has 
on controller performance. 

Lastly, the results of a simulation study, for the design of a station keeping 


controller that will be used on the NPS Phoenix AUV, will be discussed. 


B. OVERVIEW OF CLASSICAL CONTROL TECHNIQUES 


The most intuitive and oldest means of eliminating the effects of a disturbance is 
to attempt to attenuate the disturbance at the source. This often translates to corrective 
measures in the system. For example, modifying the electronics in a sensor so that the 
noise is reduced, is one common application of this technique. Other examples are 
reducing friction forces in a servo by using better bearings, or moving a Sensor to a 
position where the disturbances are smaller. Although this method of reduction at the 


source is beautiful in its simplicity, it is often impossible to achieve. 
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1 Feedback Control 


If the disturbances cannot be rejected at the source, feedback control can be used. 
For this method, the manner in which the disturbance enters the system must be known, 
and the system must be both controllable and observable. In this way, the effects of the 
disturbance can be mitigated by using local feedback. 

The classical control problem, Figure 5.1, simply stated, is, given a plant model 
Gp, design a controller, G,,, such that the closed loop system 
1. is stable and exhibits some level of robustness against plant parameter variations; 
'2. accurately tracks the reference input signal, r; and 
3. rejects the disturbance d, and noise v. 
The solution to this problem is accomplished, in general, by selecting a controller such 
that a high loop gain is obtained over the frequency range of the input disturbance, while 
obtaining a low loop gain over the range of the frequencies associated with the 


measurement noise. 


Figure 5.1 Block Diagram of a Feedback Controller 


A proper design can be obtained by simultaneously solving the two performance 


criterion, 
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over the range of frequencies of interest. This approach is typically referred to as "loop- 


(5.1) 


shaping,’’ and any of the standard control design techniques can be used to accomplish 
this goal, [Ogata 1990]. 

To eliminate any step or bias error (or to follow a step input), an integrator is 
added to the system The integrator directly changes the system's sensitivity to 
disturbances, specifically targeting bias disturbances. Integral control, in application, 
needs to be managed carefully using anti-windup methods. This compensation technique 
does not specifically reject the disturbance, but acts to alter the resulting dynamics of the 
system when a disturbance is present. 


Zs: Feedforward of Directly Measured Disturbance 


If a disturbance can be measured or estimated, feedforward control is a useful 
method of canceling its effects on the system response. Unlike feedback control, this 
method is advantageous in that it is implemented by approximately compensating for 
disturbances before they are sensed. In feedforward control, a signal from a measurable 
disturbance maybe used to generate an appropriate control force to counteract or mitigate 
the effects of the disturbance. It minimizes the magnitude of the output for the 
disturbance input without the use of error integration. 

Feedforward control alone can minimize transient errors but there are no 
guarantees of its accuracy due to its open-loop nature. Thus, feedforward alone is 
unrealistic for most applications with unsuitable open-loop dynamics. For this reason, 
feedback control is often used together with feedforward control to compensate for the 


latter's inaccuracies in minimizing the error. This approach is shown in Figure 5.2. 
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Figure 5.2 Feedforward-Feedback Controller Block Diagram 


C, MODEL BASED CONTROL 


The state-space representation of a linear system may be defined by the following 


set of equations: 
x(t) = Ax(t)+ Bu(t) 


: (52) 
y(t) = Cx(t) + Du(t) 


where A, B, C and D are determined from the differential equations describing a system. 


A block diagram depicting this system is shown in Figure 5.3. 
dt : 2 
x(t) 
se 8G : 
u(t 
i y(t) 


Figure 5.3 Block Diagram of State-Space System 


ae 


1. Linear Quadratic Regulator Control 


Performance indexes are a way of obtaining desirable output regulation without 
requiring excessive input signals. One of the most common and useful is the LOR 


performance index defined by 
oy 
1 E =x" Mx+ [(x’Qx+u" Rujic (5.3) 
0 


where M and @Q are real, symmetric positive semi-definite matrices and R is a real, 
symmetric positive definite matrix. M is called the terminal penalty matrix, @ is the state 
weighting matrix, and RK is the control weighting matrix. 
Using a performance index of this form, subject to the system dynamics given in 
Equation 5.2, a linear state feedback (LSF) control law of the form 
u(t) = Kx(t) (5.4) 
was found to be optimal for minimizing Jzor. This control law results in the block 


diagram depicted in Figure 5.4. 


+ 
a 
+ 


y(t) 





Figure 5.4 Block Diagram of Closed-Loop State-Space System with LSF 


If a unique positive definite solution to the steady-state matrix Riccati equation, 


O-PBR"B'P+PA+A'P=0, (5) 
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represented by P, existed, then the LSF control law in Equation 5.4 results in the 
minimization of Jror. The minimization of Jzgr implies a desire to minimize both 
excessive output excursions and the control effort required to prevent such excursions. 
The adjustable weights M, Q, and RK can be used to obtain an appropriate compromise 
between these two conflicting goals. The optimal control for this problem then becomes 
u(t) =—-R™"B’ Px(t). (5.6) 
Use of LQR control results in the optimal gain K and optimal pole positions. This 
method works for time-invariant or time-varying systems and is just as easy for MIMO 
systems as for SISO systems. Like the classical control methods previously discussed 
however, this control method is a feedback approach which attempts to tailor the system 
response while not specifically rejecting a disturbance. However, if the state dynamics 
matrix is augmented, to include the internal states of the disturbance signal, the 
disturbance can be rejected. 


2 LQR Control with Disturbance F eedforward 


Given the state-space equation of the system, 


x=Ax+Bu+F,d (5.7) 


if B and F are collinear, then the system can be rewritten as 
x =Ax+B(u+ae), (5.8) 
and the contro] law can be written in the form 
u=-Kx- ad, (5.9) 
as long as the disturbance is measurable. 
However, if B and F are not collinear, then direct feedforward cannot be used and 
a disturbance estimator must be employed. With this approach, the state space system 
must be augmented with a disturbance model, and a controller designed based on this 
augmented system. With a design of this form, the LQR controller can account for the 
effects of the unwanted input. This can be done only when some assumptions about the 
form of the disturbance model can be made. 
The disturbance state z, with internal dynamics Ag, may be represented by 


Z=A,Z, (5.10) 
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where 
d=C,,Z. Galt) 
Augmenting the system states, x, with the disturbance states, z, yields the following state 


space equation 


x(t A FUGA xX B 
, a = aa |] | u(t), (i2) 
z(t) 0 A, z(t) 0 
which may be expressed compactly as 
xX, Q=A,x, O+B u+w; u=—-K,x-K,z, e5) 


with K; and K2 representing the feedback and feedforward gains, respectively. 

With the disturbance dynamics included in the state dynamics matrix, the cost 
function, Jigr, is minimized to determine a new set of gains based on the augmented 
system. A disturbance state estimator is necessary to provide the internal disturbance 
states and the LQR controller uses these estimated states in the state feedback loop. 
Figure 5.5 depicts this system. Note, that the augmentation of the states has no effect on 
the dynamics of the disturbance estimator since it is uncontrollable from u(t). The 
disturbance state estimation is driven by measurements from the output or disturbance as 
available. With the formulation shown in Equation 5.13 the effects of the disturbance can 
be reduced, and in theory cancelled if a perfect disturbance model is available. 

To this point, the control methods discussed have been based on linear or 
linearized systems. The ability to use these tools for the design of controllers that will be 
used on nonlinear systems is some what limited. In general, the linear systems may not 
very robust to model mismatch which can result in system instabilities, although, 
robustness measures can be implemented into the design process by H2/H. and LMI 
techniques, [Silvestre 1998a,b]. 

3. Nonlinear Methods 


In nonlinear control, the concept of feedback plays a fundamental role in 
controller design, as it does in linear control. However, the importance of feedforward is 
much more marked than in linear control. Feedforward is used to cancel the effects of 


known disturbances and provide anticipatory corrections in tracking behavior. Very 


9D 


often, it is impossible to control a nonlinear system without feedforward compensation. 
Note that a model of the plant is always required for feedforward compensation, although 


this model need not be very accurate. 


d(t) 





eS q ce d(t)+v(t) 
ae 
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Figure 5.5 Block Diagram of Closed-Loop State-Space System with LSF and Estimated 


Disturbance Feedforward 


There is no general method for the design of nonlinear controllers. What is 
available to the designer is a collection of tools that are applicable to particular classes of 
nonlinear control problems. These nonlinear design tools can be placed in one of five 
categories: Trial-and-error, Feedback/Input Linearization, Adaptive Control, Robust or 
Sliding Mode Control and Gain-scheduling. Unlike the linear control discussion, the 
sections that follow will briefly discuss some of the salient points of these techniques. 
For further detailed information on each of these design tools, the reader is referred to 


{Slotine and Li 1991}. 
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Trial-and-error can be used to develop controllers. This method is similar in 
approach to linear lead-lag compensator design using Bode plots. The goal is to use 
analysis tools, (1.e., phase plane, describing function, Lyapunov analysis), to assist in the 
search for a control solution that can be qualified by analysis and simulation. Experience 
plays a major role in this technique, which, for complex system typically fails. 

Feedback linearization deals with techniques for transforming complex models 
into equivalent models of a simpler form, [Slotine and Li 1991]. In this nonlinear design 
methodology, the idea is to first transform the nonlinear system into a full or partially 
linear system. Once this has been done, then any of the linear design tools may be used 
to develop the necessary control system. Two draw backs of this method are that it 
typically requires full state feedback and it is not very robust to parametric uncertainty or 
disturbances. These draw backs can be overcome by the use of either robust or adaptive 
control methods. 

For uncertain or time varying systems adaptive control is very useful. Current 
adaptive control design applies mainly to systems that have well known dynamics, but 
unknown or slowly varying parameters. Adaptive controllers, whether developed for 
linear or nonlinear systems are inherently nonlinear. For nonlinear systems, adaptive 
control can be viewed as an alternative to robust nonlinear control. 

Robust nonlinear control techniques have proven very effective in a variety of 
practical control problems, [Healey 1993, Marco 1996, Yoeger 1991, Young 1996]. The 
controller is designed based on consideration of both the nominal model, and some 
characterization of the uncertainties associated with the model. Sliding Mode Control 
provides a systematic approach to the problem of maintaining stability and performance 
in the presence of modeling inprecisions. 

Gain scheduling is an attempt to apply linear control methods to the control of 
nonlinear systems. It was originally developed by the aircraft industry for the control of 
high precision aircraft. The idea is to select a number of typical operating points which 
cover the systems range of operation. The plant is then linearized and a controller 
designed for each of these points. Between operating points, the gains of the 


compensator are scheduled resulting in a global controller. The main problems with gain 


oF 


scheduling is that it has only limited stability guarantees for nonlinear operations, and the 


computational burden of computing many linear controllers. 


D. CASE STUDIES IN DISTURBANCE REJECTION 


This section will present and discuss three distinct cases of disturbance rejection 
for three different disturbance.inputs. The performance of each of the control designs 
will be evaluated by using the nonlinear EOM, Equation 2.94, in simulation studies. 
During the simulations, each control design will be subjected to a simple harmonic 
disturbance input, a PM spectrum based disturbance input and real disturbance data 
obtained from Monterey Bay. Sample data records for each of these disturbance inputs 
are shown in Figure 5.6. 

The three cases are summarized below: 

e Case I: High Gain LQR Control. Equation 2.94 will be linearized around a 


nominal operating point. Based on this linearized model the control gains for 
a full state feedback controller are calculated using a LOR method. 


e Case II: LQR Control] With Estimated Disturbance Feedforward. Employing 
the linearized model from Case I, the system is augmented with an AR model 
representing the disturbance dynamics. The augmented system is used to 
calculate the control gains, and the AR model is used as a basis for a 
disturbance estimator. This controller uses full state feedback, with estimated 
disturbance state feedforward in the control calculation. 


e Case III: Sliding Mode Control (SMC) With Measured Disturbance 
Feedforward. A model based sliding mode controller will be developed which 
embeds the disturbance in the control formulation. This controller relies on 
full state feedback with measured disturbance feedforward. 


1. Case I. High Gain LQR Control 


Using the | DOF surge EOM as a model, it is necessary to linearize this system of 
equations in order to use linear techniques in the controller design. Linearization can be 
performed a with a variety of approaches, stochastic [Leira 1987], harmonic 
[Heyns 1995] or nominal operating, pointwise linearization, condition [Riedel 1998a], 


but for this design, since the control law will attempt to allow an AUV to hold position 
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Figure 5.6 Sample Disturbance Input Time Series 


with zero ground velocity, the system is linearized about the steady state solution to 
Equation 2.94 while in the presence of a steady current. Performing this steady state 


analysis yields, 








le ose = au, ? (3.14) 
2 
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as the nominal operating point, where n, must be real and the same sign as u,,. Using a 
standard Taylor series linearization, the linearized system of equations, in state space 


from, can then be written as 








x 0 ] Ome 0 ] 
u,|=|O 2a@sgn(u,,)u,, 1 jiu, i+ 0 n+|O|u, 
F 0 if i. SLE ue u,,sign(n, )+ en )n, 0 
q Ta co eo 
I O O|j x 
y=|0 1 Oflu, 
OO el a 


Ons 

If it is assumed that the vehicle will be operating in a -0.1 m/s steady current, and the 
parameters @ f, yand Tare available then Equation 5.15 can be evaluated numerically. 

Since the parameters identified in Chapter IV were obtained from a discrete filter, 

and it is desired to implement the to be developed control law in a digital computer, the 

state space equations must be converted into a discrete form. Using a standard Euler 


discretization, Equation 5.15 can be represented in discrete form as 


x(k +1) 0 dt 0 x(k) 
u(k+1)}/=|O0 1-2a@sgn(u,,)u, dt dt TON) | toes 
FR*V] |o Lin tae 1 hae JL 
L Tt Td 
0 dt 
0 n(k) +} 0 Ju,(k). (5.16) 
(Fu,, sign(n, )+ 2 sen(n, )n, at 0 
Il O O}j x(k) 
y(k)=|0 1 Off u,(k) 
O O 11 Fik) 


Using standard optimal control techniques the solution for the optimal (LQR) 


controller can be found as 
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n=-R'B' Sx (5.17) 

where S is found by solving the steady state algebraic Riccati equation (ARE) 
A’S+SA-—SBR’B'S+Q=0, (5.18) 
for the positive definite matrix S. In Equation 5.26, @ is the weighting matrix on the state 
error, and R is the weighting scalar, since this is a single input system, that invokes a 
penalty against the contro] effort. The LQR approach will always yield a stable system, 
as long as the Riccati equation provides a positive definite solution matrix S$, for which 

the system must be controllable and full state feedback available. 

The following sections will show the simulated performance of the LOR 
controller when subjected to various disturbance inputs. The purpose of the simulations 
in these sub-cases is to provide a baseline by which to compare the performance the 


controllers developed in Cases II and III. 
a) Monochromatic Disturbance Input (Case Ia) 


Using the sine wave disturbance input depicted in Figure 5.6, the LQR 
controller was formulated and simulated for various control weighting values (R). A plot 
of the position response for one of these simulations is shown is Figure 5.7. This 
response, the "best" of the many simulations, is the result of a high gain controller. As 
can be seen, the oscillations about the commanded position (O meters) are significant and 
poor disturbance rejection is obtained. In addition, there is an obvious offset caused by 
the mean disturbance. This offset can be corrected by incorporating integral control into 
the LQR design, however, integral control will not correct the severe positional 
oscillations. 

During these simulations, the propellers were not limited, i.e., no 
Saturation. The control input necessary to obtain the positioning shown in Figure 5.7 is 
displayed in Figure 5.8. As shown is this figure, the propeller oscillations are extreme 
considering that the model and parameters used in the simulations are based on the NPS 
Phoenix AUV which has a maximum propeller revolution of 800 rpm. Once again, it 
must be pointed out that the purpose of the studies in Case I is to obtain a baseline for 


comparison. 
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Figure 5.7 Position Response for Case Ia with Monochromatic Disturbance Input 
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Figure 5.8 Propeller Response for Case Ia with Monochromatic Disturbance Input 
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b) PM Spectrum Based Disturbance Input 


Using the same three state model as with Case Ia, simulation studies using 
a PM spectrum based disturbance input were conducted. The input disturbance was 
based on a significant wave height of 1 meter in a water depth of 45 meters. The vehicle 
was assumed to be operating at a 25 meter depth. The goal of this simulation study was 
to determine the control performance based on a disturbance input which contained a 
range of frequencies which the vehicle may encounter. 

Using the controller design resulting in the responses displayed in Figures 
5.7 and 5.8, a simulation was conducted resulting in the position response shown in 
Figure 5.9. In this particular simulation, the standard deviation of the position response is 
significantly reduced due to the magnitude of the disturbance input. Comparing 


disturbance inputs between Cases Ia and Ib, it can be seen that the magnitude of the PM 
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Figure 5.9 Position Response for Case Ib with PM Spectrum Based Disturbance Input 


based disturbance is 15 times less than that of the monochromatic disturbance. This 


reduction in oscillation magnitude is also reflected when comparing position responses. 
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The propeller input response shown in Figure 5.10, is also significantly less, but still 


exceeds the maximum propulsion system input of 800 rpm. 
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Figure 5.10 Propeller Response for Case Ib with PM Based Disturbance Input 


Cc) Monterey Bay Disturbance Input 


Using the control design that resulted in the responses displayed in Figures 
5.7-5.10, a third set of simulations was conducted. In this set of simulations, transformed 
wave buoy data obtained in Monterey Bay, CA was used as the disturbance input. This 
data was obtained from a Datawell® Waverider Buoy deployed April 9, 1998, from the 
research vessel R/V POINT SUR, during an NPS oceanography class (OC4610) cruise, 
under the direction of Prof. Thomas Herbers. 

The wave buoy, according to Defense Mapping Agency navigation charts, 
was deployed in approximately 45 meters of water. The wave elevation data obtained 
from the buoy was transformed to a subsurface velocity record, at a depth of 25 meters, 
using the procedure outlined in Chapter IV. A sample of the resulting time series, with a 


-0.1 m/s steady current superimposed, was displayed in Figure 5.6 
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The position response using this input disturbance is shown in Figure 5.11. 
The standard deviation of this response is approximately one-half of the response 
obtained from the PM based disturbance input with the same control design. Once again, 
this is due to the fact that the standard deviation of the Monterey Bay input disturbance is 


about one-half the PM disturbance. 
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Figure 5.11 Position Response for Case Ic with Monterey Bay Disturbance Input 


It is interesting to note, that although the position response reduced by a 
factor of two, when compared to the PM based case, the propeller input response did not, 
see Figure 5.12. This is due to the fact that the frequency content of the input 
disturbances is much different. This is evident by referring back to Figure 5.6. The Bay 
data contains more high frequency components causing the propulsion system to respond 
much more. 

Since the propulsion system response is still in excess of maximum output, 
tuning of the controller gains must be performed to bring the propeller rpms within limits. 
By adjusting the input weighting scalar R, and reducing the controller gains the maximum 
commanded propeller revolutions can be reduced as well as reducing the sensitivity of 


the controller to high frequency "noise.’’ This reduction of propeller input is at the 
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expense of increased position error. These results are shown in Figures 5.13 and 5.14. 


As shown in Figure 5.13, the standard deviation has increased by a factor of two in order 


to keep the propeller revolutions within propulsion system limits. 
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Figure 5.12 Propeller Response for Case Ic, Monterey Bay Disturbance Input 
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Figure 5.13 Position Response for Case Ic, Monterey Bay Disturbance Input 
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Figure 5.14 Propeller Response for Case Ic, Monterey Bay Disturbance Input, rpms 
Within Design Limits 

Figure 5.15 shows graphically the relationship between the level of control 
input and the level of disturbance rejection for the standard LQR solution subjected to 
real wave data from Monterey Bay. The position covariance is normalized by the 
covariance of the "free floating" or uncontrolled response of the vehicle, and the input 
covariance is normalized by the maximum rpm available from the propellers. This 
analysis can give a "feel" for how tight a control law must be provided to achieve a 
reasonable disturbance rejection. 


2 Case IT. LQR Control with Disturbance Estimation Feedforward 


The problem that now must be addressed is how to achieve better performance. It 
has been shown , that by embedding an estimator of the disturbance into the control 
system design, improved performance may be obtained [Grimble 1995, Riedel 1998a]. 

| As outlined in Chapter If, an AR model of the wave disturbance may be written 
in state space form as 
X (kK +I =A,X (kK) + B vk) 


Se, 
yy (kK) =u, (A) HC, xX, (Kk) oe 
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Figure 5.15 Comparison Of Control Input Covariance To Normalized Vehicle Position 


Covariance, Monterey Bay Wave Data 


Augmenting the vehicle state equations with the disturbance state equations, a 
new control law may be developed using the estimated disturbance states. Defining the 
new state vector as 

XY = i oo (5.20) 
aug xX. 
where the disturbance states are given as 
X,=[X (kK+N-)),....X, (K+) Ff , (52 
the new control law may be designed, using the separation principle, assuming all states 
are measurable. As in the previous optimal control discussion (Case I), the ARE is 
solved to obtain the appropriate gains for the augmented system. This augmented system 
is represented by 


0 A 


w 


A ehicte ia ce B vehicle 0 
X wae (K +1) = X aug (R)+)“O In(k)+] |v (kD. (5.22) 


108 


With the control law determined, the estimator must be designed. Using optimal 


estimation theory, an estimator of the form 
X (k+l) =(A, -LC, )X,, (k) + Lu, (k), (5.23) 


where uk) is the current disturbance measurement, is developed. This estimator is used 
in conjunction with the control law developed, and its implementation, in block diagram 


form, is represented by Figure 5.5. 
a) Monochromatic Disturbance Input 


To display how well this design procedure can work if an accurate model 
of the disturbance is available, consider the case of the monochromatic input disturbance. 
Since the precise model of this disturbance is known, when this model is embedded in the 
control system design, perfect cancellation of the wave disturbance effects on vehicle 
positioning may be obtained. These results are displayed in Figures 5.16 and 5.17. Now 
these results are for demonstration purposes only, and perfect cancellation of the wave 


disturbances will not be possible since an exact model of the random sea is not available. 
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Figure 5.16 Position Response for Case Ila, Monochromatic Disturbance Input 
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Figure 5.17 Propeller Response Case Ila, Monochromatic Disturbance Input 


Although the propeller response is almost identical to the results displayed 
in Case Ia (Figure 5.8), by having an estimate of the disturbance states to feedforward the 


propeller input is properly phased to cancel the disturbance. 


b) PM Spectrum Based Disturbance Input 


Using this approach of a model based disturbance estimator with a LOR 
controller appears to be an excellent method of canceling the disturbances acting on an 
underwater vehicle, that 1s if the model of the disturbance is known. If the exact model 
of the disturbance is not known, the question is; Is improved disturbance rejection with 
this method possible? 

Adopting the AR modeling techniques presented in Chapter IV, a sixth 
order AR model for the PM based disturbance was developed. Using this linear model of 
the disturbance dynamics and the same input weighting scalar as was used in Case Ib, a 
combined controller/estimator was developed. Using this developed compensator in 
simulation, improved performance was observed, see Figure 5.18. The position response 


with the augmented disturbance model improves by a factor of 1.5. 
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Figure 5.18 Position Response for Case IIb, PM Based Disturbance Input 


The improvement in propeller input follows the trend displayed in the 
comparison between Cases Ia and Ila. The standard deviation of the input response has 
not changed significantly, as evident in comparisons between Figures 5.19 and 5.10. 
What has changed is the control input phasing, again due to the disturbance feedforward, 
thus allowing this design method, even with a low order disturbance model, to obtain 


improved disturbance rejection. 
Cc) Monterey Bay Disturbance Input 


Using identical weighting values that went into the design of the control 
laws used in the simulations presented in Case Ic, Figures 5.11-5.14, and a sixth-order 
AR model representing the Monterey Bay disturbance, improved performance was again 
realized. As can be seen in Figure 5.20, there is a 150% improvement in station keeping 
as compared Case Ic, and the control input requirements are significantly less, see Figure 
5.21. Although the standard deviation of the commanded control input is well within the 
maximum revolutions able to be provided by the propulsion system, there are some 
inputs which exceed the limit of 800 rpm. In order to bring the commanded control input 
within limits, as before, the control gain must be reduced. 
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Figure 5.19 Propeller Response for Case IIb, PM Based Disturbance Input 
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Figure 5.20 Position Response for Case IIc, Monterey Bay Disturbance Input 


With the control gains adjusted so that the commanded control input remained 


within propulsion system limits, the positional error increased by a factor of two, while 


EZ 


the control input reduced by a factor of three The results of this tuning are shown in 


Figures 5.22 and 5.23. This reduction in control effort is particularly important given the 


fact that power consumption is the downfall of AUVs. 
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Figure 5.21 Propeller Response for Case IIc, Monterey Bay Disturbance Input 
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Figure 5.22 Position Response for Case IIc, Monterey Bay Disturbance Input, rpms 
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Figure 5.23 Propeller Response for Case Ic, Monterey Bay Disturbance Input, rpms 
Within Design Limits 


3. Case III. Sliding Mode Control with Measured Disturbance 
Feedforward 


Beginning with Equation 2.94, a sliding mode controller was formulated using 
standard SMC techniques, [Slotine 1991]. The sliding surface o was defined as a 


function of the position error, 
o-(4 +A) (ema ), (5.24) 
dt a 


and the time derivative of o was defined as 

Go =-nsat(o/¢). (5.25) 
By defining the sliding surface in this manner, stability is guaranteed, based on Lyapunov 
analysis, since 

oo <0, Vt 0. (5.26) 
Taking the time derivative of Equation 5.24 and equating it to Equation 5.25, the control 


input may be determined. 
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Using the signed square root of Equation 5.27, the commanded control input is found. A 
detailed description of this controller design approach may be found in [Riedel] 1998b]. 

| As seen in Equation 5.27, the commanded control input is a function of the 
. system states, the fluid velocity (including the first and second derivative), the command 
inputs, and the "to-be computed" control input, due‘to the fact the system represented by 
Equation 2.94 is non-affine. To compute the required control input requires solving a 
difference equation in 1, as well as measurements of the fluid velocity and its first and 
second derivative, making this control law extremely complex and possibly difficult to 


implement in real-time. To overcome these difficulties, some simplifications need to be 


made. 
a) Monochromatic Disturbance Input 


If the thrust reduction term is ignored and treated as an unmodeled disturbance, 
Equation 5.27 reduces to, 


+ Fsign(u,)+... 








6 —2au, lou. u, 


eye ae (5.28) 








(ale Oey 
nhl == me! eves ss —2A(Qu, u, 


U ~X om -A (lu, oie) 

which requires only system states, fluid disturbance measurements and command inputs. 
To display how well this controller is capable of performing, again, consider the case of a 
monochromatic sine wave disturbance input, where the disturbance and its first and 
second derivative are known. When direct feedforward of the measured wave disturbance 
is embedded in the control system design, perfect cancellation of the wave disturbance 
effects on station keeping may be obtained. The simulated response of the PHOENIX, 


initially at five meters and closing to a commanded range of 0.5 meters, is displayed in 
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Figure 5.24. Again, the results in Figure 5.24 are for demonstration purposes only, as a 
comparison with Case Ila. Perfect cancellation of the wave disturbances is not expected, 


since exact measurement of each wave disturbance component, 1.e., u eo u ¢ and u pe is not 


possible. The interesting point in this case is the propeller response. 


6 





Figure 5.24 Disturbance Cancellation Case IIIa, top to bottom respectively, position vs. 


time, propeller RPM vs. time, and a phase plane plot of the sliding surface 


When comparing the propeller response between the three cases that used a sine 
wave disturbance input, it can be seen that the SMC (Case IIIa) by far out performs the 
other designs. The position response is as desired, perfect cancellation, and the 
propulsion system is well within limits. This result is due to the fact that the system 


attempting to be controlled is highly nonlinear, requiring a nonlinear controller. 


b) PM Spectrum Based Disturbance Input 


Prior to continuing with any simulations to determine positioning 
performance, several simulations were conducted to determine the performance that 
could be obtained from the controller with and without all disturbance components 


available. Since the PM based disturbance input was generated using the techniques in 
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Chapter IV, the derivatives and phasing were known. Based on this, comparative 
simulations were conducted between controllers which used all the disturbance states and 
ones which used only the measurable fluid velocity state for disturbance rejection. 


Comparative results of one simulation are shown in Figure 5.25. 
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Figure 5.25 Controller Performance Comparison, For A Controller That Uses All The 
Disturbance Components (Dashed Line), And A Controller That Uses Only Fluid 


Velocity For Disturbance Cancellation (Solid Line) 


As can be seen in Figure 5.25, the station-keeping improvements associated with 
including all components as opposed to including only the fluid velocity component is 
very small. In each case, the propulsion system response was within the vehicle's 
capability. As a result of the comparisons, it was determined that by using only the fluid 
velocity measurements, significant improvement with regard to positioning may be 


achieved. Based on this, the SMC takes the form 
+F |sign(u,)+... 








O-2Qu. lou. u, 





n|n| = — Be 2 A(aa le ee . (5.29) 
B\t 3 


y-A (u, +1 — eee) 


X com 
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which will be used for all remaining simulations. 

Using the PM based disturbance input allowed the SMC to be tuned to so 
that the controller would meet bandwidth requirements, limit propulsion system 
oscillations and avoid chattering. Controller parameters which provided a balanced 
design consisted of 7 = 100, A = 1.0, and ¢= 0.5. The simulated position response of the 


Phoenix conducted with this SMC design is shown in Figure 5.26. 
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Figure 5.26 Position Response for Case IIIb, PM Based Disturbance Input 


The position response shown in Figure 5.26 has a standard deviation of 
6.4 cm. This is twice as much as Case IIb, Figure 5.18, however, the standard deviation 
of the propeller input for the SMC design is one-half a large as the LQR with disturbance 
estimator design, and is always within propulsion system limits. This can be seen in 
Figure 5.27. In addition, when comparing the two propeller responses, it appears that the 


SMC has a smoother output which will extend the life of the propulsion system. 
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Cc) Monterey Bay Disturbance Input 


Using the same design parameters that allowed the controller to achieve 
the performance depicted in Figure 5.27, the system was simulated with the disturbance 
input obtained from Monterey Bay. Since the disturbance magnitude of the Monterey 
Bay data is less than the PM based input, it 1s expected that the position response would 
also be less. By referring to Figure 5.28, it can be seen that this is in fact the case. 

The standard deviation of the position response has improved over the 
LQR based controller (Case IIc) by a factor of 1.4 with only a slight increase in propeller 
rpms (5% compared to a maximum of 800). These results are shown in Figure 5.29. 


4. Disturbance Rejection Case Comparison 


After conducting the simulations for each of the cases with the various 
disturbance inputs, it was apparent that Case III, the SMC with measured disturbance 


feedforward, out performed the other two cases and to most this is no surprise. What is 


interesting is the amount by which it out performed the other cases. 
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Figure 5.27 Propeller Response for Case IIIb, PM Based Disturbance Input 
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Figure 5.28 Position Response for Case IIIc, Monterey Bay Disturbance Input 
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Figure 5.29 Propeller Response for Case IIIc, Monterey Bay Disturbance Input 


Using the PM based disturbance input, simulations were conducted for each of the 


three case previously analyzed, LOR, LQR with disturbance estimator, and SMC with 


120 


disturbance measurement feedforward, with gains ranging from high to low. The attempt 
was to reproduce the "optimality" curve, Figure 5.15, for each controller to study the 
performance of each control solution. 

Conducting this study led to some very interesting results which are shown in 
Figure 5.30. As seen in this plot, the curves for each controller do not have the traditional 
optimal curve shape, i.e., aS contro] input increases position error decreases. In fact, the 
curves indicate that for control designs ranging from low to medium gains, regardless of 


controller type, it would be better to have no control at all. 
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Figure 5.30 Comparison Of Controllers For Various Gains, PM Based Disturbance Input 


The explanation for this can be seen in Figure 5.31. which superimposes 
the closed-loop vehicle frequency response, disturbance input to vehicle position output, 
over the disturbance spectrum for three different control gains, namely low medium and 
high. In Figures 5.30 and 5.31, point/curve "1" corresponds to a low gain, 
point/curve "2" to a medium gain solution and point/curve "3" is high gain control. 

As Figure 5.31 displays, the low and medium gain solutions, with this 
particular disturbance, actually excites the vehicle, and not until a high gain solution is 


implemented does the vehicle actually reject the disturbance. Using this as an analysis 
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tool, the range of acceptable gains, for a particular disturbance input may be determined. 
In addition, it was quite evident that by feeding forward the measured disturbance using 
the SMC solution, that significant disturbance rejection was capable with much less 


power consumption. 
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Figure 5.31 Vehicle Frequency Response, (Disturbance Input To Position Output), 
Superimposed Over The PM Based Disturbance Input 


EK. SUMMARY 


This chapter has outlined the various disturbance rejection techniques available to 
the control engineer. It has provided a summary of classical, modern and nonlinear 
control methodologies. Three case studies, which represent the basic design methods 
used to reject disturbances, were conducted and discussed for three different disturbance 
inputs. These studies showed that the SMC with measured disturbance feedforward is a 


far superior design approach for this particular class of problem, and provides significant 


2 


disturbance rejection performance for the same input power. Finally, an analysis 
approach that may be used to study gain selection and performance estimates has been 


presented. 
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VI. DISTURBANCE COMPENSATION CONTROLLER (DCC) 


A. INTRODUCTION 


This chapter will discuss the development of the real-time disturbance 
compensation controller (DCC) which will allow an AUV to dynamically position itself 
in the presence of waves. The chapter will begin with an overview of the DCC, followed 
by a discussion of an asynchronous Extended Kalman Filter for state and disturbance 
estimation. This nonlinear estimator is critical to the DCC performance since the SMC 
requires full state feedback, and not all states are measurable. In addition, the EKF 
provides the controller with a smoothed estimate of the unmeasured fluid velocity which 
is used to compensate for the wave induced disturbance. 

Next, through the design and implementation of an asynchronous simulator, 
which realistically models the vehicle dynamics, the sensors including noise and the 
sensor processes, the DCC is tuned and the achievable performance is demonstrated. 

Lastly, it is shown that by properly weighting the noise covariance in the 
estimator the DCC reduces the transmission of sensor noise into the propulsion system 


while still maintaining the ability of the vehicle to hold position. 


B. DCC OVERVIEW 


The design of the disturbance compensation controller can be looked at as an 
optimization problem since there are competing goals. First, since the design 
requirement is to minimize position error in the presence of disturbances, a high gain 
control, as Chapter V discussed, is desirable. Using high gain control, the system 
becomes sensitive to measurement noise and uncertainty, thereby requiring the gain to be 
reduced to maintain stability. 

The estimator is needed to provide the unmeasurable states to the controller, and 
to filter the sensor noise thereby improving the systems performance. Here, the 
requirement is to accurately track the signal, again requiring a high filter gain, while 


smoothing the noise, (a low gain). As with the controller, trade-offs must be made. 


ZS 


The overall goal is to develop a combined controller/estimator which, when 
implemented, will enable the vehicle to maintain position while using noisy sensor 
information. The output of this system is a commanded voltage that is sent from the 
DCC process to the real-time execution computer, without excessive lags to ensure 
stability. A mathematical description to the above problem is given below, with a block 


diagram of the DCC in provided in Figure 6.1. 


State : x =|X.a P|, d=u, 
System : x= f(x,n,d); y’ =|x,u,,u,,|=Cx+Dd 
Disturbance: X,=Ax,+V uu, =Cx, ; (6.1) 


an 


Estimator : Ika, |= EKF( f(i,n,d),.A.y.n) 


Control law: n= smc(x,U,,x 
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Figure 6.1 Block Diagram of Disturbance Compensation Controller (DCC) 


c. STATE AND DISTURBANCE ESTIMATION 


There are many methods available to estimate states and disturbances in practice 


today. A few of these include the Luenberger Observer [Ogata 1990] and the Kalman 
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Filter [Gelb 1974] for linear systems, and the Sliding Mode Observer [Canudas De Wit 
1991], the "Rajamani" Observer [Rajamani 1998] and the Extended Kalman Filter for 
nonlinear systems. Each method has both pros and cons depending on the application. 
For this work, an Extended Kalman Filter was chosen since a relatively accurate vehicle 
model is available, and the stochastic nature of the disturbance. 

Kalman filtering is the process of recursively updating an estimate of systems 
states based upon measurements corrupted by noise. The system state is a collection of 
variables that describe the dynamics of a system, and in this case they are position, 
relative velocity and propeller thrust, of which only relative velocity is measurable. 

System states are updated with knowledge of system dynamics (vehicle model), 
measurement dynamics (measurement model), system noise (modeling uncertainty) and 
measurement noise (measurement errors). The system model is not perfect in describing 
the dynamics of the vehicle and will contain a certain amount of uncertainty, called 
system noise. There is also some uncertainty associated with each measurement taken. 
- This uncertainty can be composed of both random white noise and a bias. Measurements 
which cannot be directly obtained, such as fluid velocity, are related to measurements 
which are directly obtainable, such as relative velocity and ground velocity, in the 
measurement model. Recursively updating means the Kalman filter does not need to 
keep record of all past measurements, only the most recent ones. 


I Model and Filter Development 


Using the three state surge model developed in Chapter II, and a four state AR 
model for the wave dynamics, an augment state and disturbance model was formed, and 
used as the basis of an EKF. This model allows the disturbance to be treated as an 
additional state, where the vehicle states and disturbance estimates are filter outputs. The 


augmented vehicle and disturbance model is given by, 


ey 


X=uU, tu, 











u.=Qu,ju.j+F 
= 
ae ino 
T T T 
MNS =Us : (6.2) 


nae =U, =X,3 

Xy3 =A,X,, a Xe 
Xy4=V 

y=[x,u,,u,]' 


where the AR coefficients are found using the procedure outline in Chapter V. 


Ze Kalman Filter Algorithm 


Using standard design techniques [Gelb 1974], the filter was developed and 
implemented using the following algorithm. First, the system model matrix A, system 
noise matrix Q, measurement matrix C, measurement noise matrix R, and the error 
covariance matrix P are initialized to appropriate values. The error covariance matrix can 
be thought of as a level of uncertainty in the state vector. Then the state vector, error 
covariance and measurement vector are propagated one time step using the model. 

When the new measurement is received, the innovation is calculated based on the 
difference between the measured values and the estimated values. Using the propagated 
error cOvariance, measurement noise matrix and measurement matrix, a gain is 
determined for the state vector and error covariance update. This process of propagating 
and updating is repeated through out the length of the vehicle mission. This recursive 


algorithm, in discrete form is given by, 


Of (x ,,.>1) 
®, 4... =exp(AT); A=———-— 
OX aug = 
X eset = On aXiine 
Pein = D/P eee ae +Q4 ° (6.2) 
G, a Ph (bebe hy, as RJ" 
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wk =X t+ G, Ly, — bX p41] 
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where @® represents the linearized system dynamics matrix, and h=C since the 
measurements are linear in the state. The continuous linearized matrices for this 


particular design are given as, 





0 l 0 0 1 0 0 
0 2a%, signté,) 1 0 0 0 0 
0 =In,,n| as 0 0 O 
U: 
i () 0 0 0 1 0 O 
0 0 0 0 0 1 =O 6.3) 
0 0 0. Sie 6. 4 
0 0 0 0 0 0 0 
1000000 
c=10 10000 0 
0000100 


3. Asynchronous Data Processing 


In the preceding discussion, the data contained in the measurements was assumed 
to be received at the same time with equal intervals through out the mission. In reality, 
all measurements are not received at the same rate, therefore, the EKF design must allow 
for this asynchronous sampling rate. In the Phoenix AUV, the vehicle control loop 
currently runs at 8 Hz, while the RDI DVL runs at 2 Hz, and the SonTek ADV at 6 Hz. 
(See Appendix B for a more through description of sensor operations). The main data 
acquisition process samples the sensor processes at the same frequency as the control 
loop, however, if the sensor has not yet updated, the data acquisition process records the 
value of the previous time step. The filter allows for the varying measurement rates by 
using a dynamic switching of the measurement matrix, C, [Healey 1998]. The 
measurement matrix basically uses a zero-order hold on the measurement channel that 


has not been updated, and propagates the state using the previous measurement. 


D. ASYNCHRONOUS SIMULATOR DEVELOPMENT 


Using the filter design from the previous section, and the sliding mode controller 


developed in Chapter V, Equation 5.37, an asynchronous simulator was developed for 


IZ? 


design validation. The simulator contains the non-linear vehicle dynamics, Equation 
2.94, asynchronous sensor models with measurement noise, seaway dynamics and the 
DCC. Using this simulator as a design tool allowed the DCC’'s control and estimation 
parameters to be adjusted prior to real-time implementation. Figure 6.2 shows a sample 
of the sensor outputs, during one of the simulation runs. As seen, the position output, 
which is a product of a navigation filter, is at 8 Hz. The relative velocity, u,, which is 
measured by the ADV, is at 6 Hz, and the ground velocity, u,, from the RDI is recorded 
at 2 Hz. In addition, the ADV output has noise imposed on the signal representative of 


the vendor advertised levels. 
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Figure 6.2 Asynchronous simulation with realistic noise models - Disturbance from PM 


Spectrum, H,= 1 meter, operating depth 10 meters 


Using this developed simulator, the DCC was adjusted to achieve an optimum 
design. The gains in both the controller and filter were adjusted so that performance 
requirements discussed earlier were met. Sample results showing the performance of 
final design are given in Figures 6.3 and 6.4. As can be seen, the estimates of both 
position and thrust track the actual values, and the position response is maintained within 


a standard deviation of 8 cm. This performance is extremely good, recalling that the 
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same controller, with the same disturbance input was able to achieve a standard deviation 


in position of 6.5 cm, without noise and using full state feedback, see Figure 5.31. 
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Figure 6.3 Simulated and Estimated Position Response, Using Final DCC Design 
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Figure 6.4 Simulated and Estimated Thrust Response, Using Final DCC Design 
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The stability performance of the estimator is shown through simulation, see 
Figure 6.5, since there are no formal proofs to determine the stability of combined 
nonlinear estimators and controllers. As seen, the error covariance levels all converge 
indicating a stable nonlinear filter design. Some of the covariance levels may appear to 
be "too high" giving the feeling that the filter is not properly designed, however, design 
decisions must be made to ensure that the filter lags are not too excessive, and that the 


estimator tracks well. 





Figure 6.5 DCC Error Covariance Evolution 


E. INITIAL IN-WATER TESTING 


Using short missions, that DCC was adjusted to achieve acceptable performance. 
These runs were performed on March 25, 1999, in Monterey Harbor. Of concern, was 
the amount of noise that was resident on the ADV sensor. This noise was far beyond the 
level which the vendor advertised. Using the design results from the simulations, the 
DCC was implemented in the Phoenix AUV. Figures 6.6-6.10 display initial results. As 


seen, the filter tracks the signals extremely well, including the noise. 
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Figure 6.6 Short Segment In-Water Results, Position for Rapv=10 
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Figure 6.7 Short Segment In-Water Results, Relative Velocity for Rapy=10 
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Figure 6.10 Short Segment In-Water Results, Propeller RPMs for Rapy=10 


This tracking of the noise has significant detrimental effects to the propulsion 
system as seen in Figure 6.10. The noise had been transmitted into the controller 
resulting in severe oscillation in the propeller response. These oscillations eventually 
lead to mechanical failure of the propulsion system shafting due to the shearing of 
connecting pins. 

Using the information obtain during this set of runs allowed the filter gains to be 
adjusted to eliminate the transmission of sensor noise into the controller. Using linear 
design techniques, the combined controller filter transfer function from ADV input to 
propeller output was formed. By adjusting the level of the measurement noise 
parameters, attenuation of the noise into the control system was accomplished. These 
results are shown in Figures 6.11 and 6.12. As the Figures shown, for the smaller noise 
estimate (Rapvy=10), the noise transmission is amplified over most of the range of the 
controller, while for the higher noise estimate (Rapv=100), the input to the controller in 


attenuated. This improvement in frequency response will reduce the propeller 
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oscillations, 


system. 


Figure 6. 


thereby minimizing the chance of mechanical failure of the propulsion 
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Figure 6.12 DCC Frequency Response, ADV Input to Propeller Output , Rapv=100 
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Using the new design values, the DCC was again tested in Monterey Harbor. The 


results of this testing are shown in Figures 6.13-6.17. 
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Figure 6.13 Short Segment In-Water Results, Position for Rapy=100 
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Figure 6.14 Short Segment In-Water Results, Relative Velocity for Rapv=100 
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Figure 6.15 Short Segment In-Water Results, Fluid Velocity Estimate for Rapv=100 
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Figure 6.16 Short Segment In-Water Results, Thrust Estimate for Rapv=100 
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Figure 6.17 Short Segment In-Water Results, Propeller RPMs for Rapyv=100 


As the result of the tuning of the DCC, the performance has improved 
dramatically. As before, the DCC maintains position extremely well, with a much 
reduced propeller response. Comparing the magnitude of the estimated fluid velocities 
between the two designs, Figures 6.8 and 6.15, it can be seen that for the same magnitude 
of input disturbance, position response has remained unchanged, but propeller response 


has reduced increasing the life of the propulsion system. 


F. SUMMARY 


The design and implementation of a new Disturbance Compensation Controller 
(DCC) has been presented. The results indicate that by using a properly tuned system, 
the ability of a tetherless underwater vehicle to dynamically position itself with minimal 
input is possible. Although no formal proof for stability is available, asynchronous 
simulations have demonstrated that the DCC is stable and provides acceptable tracking 


and estimation of state and disturbance inputs. 
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VII. ESTIMATION OF WAVE DIRECTIONALITY FROM A MOVING 
PLATFORM 


x INTRODUCTION 


This chapter will outline the underlying principles used in identification of wave 
directions from standard wave following buoys. It will present the mathematical 
formulas used in determining the wave direction as a function of frequency. Extension of 
these algorithms to subsurface velocity sensors will be made, where, through the use of 
the Doppler equation, a moving AUV can be used to determine wave directions. Lastly, 
it will be shown how a control command can be obtained from the frequency dependent 
wave direction estimates. 

The information in this chapter is not new, only the application to which this 
method is applied. For more detailed descriptions of the mathematical formulations 
presented in this chapter, the reader is referred to [O'Reilly 1996} and the references 


therein. 


B. WAVE SPECTRA AND DIRECTIONAL ESTIMATES 


As discussed in Chapter II, a wave record, 7(t), measured at a fixed location can 
be represented as a linear superposition of waves of different frequencies. Wave 
components with different frequencies are usually assumed to be statistically independent 
because they are generated by random wind forces at different locations. From the 
central limit theorem it follows that the probability distribution of 7(t) is approximately 
Gaussian, consistent with many observations, [Soong 1993]. 

The procedure presently employed by many of the operational data buoys in 
based on Fourier analysis. In Fourier analysis it is convenient to work with complex 
exponentials rather than sine and cosine functions, there fore using the relation 


exp(i(@t+@)) +exp(—i(at + @)) 


cos(a@t +) = 5 (7.1) 
the expression for the surface elevation can be written as 
n(t)= > A, exp(iar), (2) 
@ 
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where 
1 
ae aoe exp(1@,, ) (7.3) 


and the summation is over both positive and negative frequencies. 


As discussed in Chapter II, the energy spectrum E(q@), is defined as 


2 
E(@)= vat) (7.4) 


where { ) indicates an average over many data records and Aq@ is the spacing of the 


frequency bands. The spectrum is closely related to the energy of the waves, and 
represents the distribution of wave energy as a function of frequency. 

To describe the spatial and temporal characteristics of the sea surface linear 
superposition of wave components is used. In exponential terms this can be represented 


as 


M(x, y,t)= > >°A,, exp(i(k(xcos@ + y sin 8) — wt)) (7.5) 
o 6 


where x, y are the horizontal spatial coordinates, and wand k obey the dispersion relation. 


The frequency directional wave spectrum is defined as 


2 
E(@) = (eal (7.6) 


and describes the distribution of energy as a function of frequency and direction. 


Cc WAVE BUOYS 


The most commonly used instrument for measuring waves in deep water is the 
"heave, pitch and roll buoy” that measures the surface height and slope in two orthogonal 
directions. The newer Datawell® Directional Waverider measures 3-component 
accelerations of the buoy which are integrated to yield the horizontal and vertical 
displacements of the buoy. The hull and mooring were designed to give the buoy good 
wave following characteristics, thereby allowing the buoy displacements to approximate 


the displacements of an actual water particle at the sea surface. 
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For a wave train propagating in the positive x-direction, the fluid particle motion 
is given by Equation 2.49. For the more general case of a wave train propagating at some 
angle relative to the x-axis, it can be shown that the flow field is given by 

u(x, y,t) =amwcos O exp(kZ) cos(k(xcos@ + y sin 8) — ar) 

v(x, y,t) =aq@sin 6 exp(kZ) cos(k(xcos@+ ysin @)-— ar) (7.7) 

w(x, y,t) = awexp(kZ) sin(k(xcos@+ ysin 8) — a) 
Let the average position of the buoy be given by x=y=z=0. For small amplitude waves, 
the expected buoy displacements are small compared to the surface wavelength, therefore 
the buoy motion can be approximated by the fluid velocity at x=y=z=0. 

For a full spectrum of waves, the buoy displacements can be expressed using 
complex notation as 

X(t)=)>) >| -iA, , cosO exp(iwt) 


o @ 
Y(t)=>)>)-iA,,, sin 8 exp(ivt) (7.8) 
@ @ 
Z(t)= >) >) - Age exp(int) 
o @ 
where the -i is due to the 90. phase difference between the vertical and horizontal 
displacements. The expressions in Equation 7.8 can be written using Fourier transforms 


as 
X(t) = >) X (@) exp(iwt) 


ri) > Y(@) exp(iwt) , (7.9) 
Z(t) = >) Z(@) exp(iwt) 


where the Fourier transforms are given by 
X(@) = >-iA,,, cos 
@ 
Y(@)=)>-iA,, sin@ . (7.10) 
@ 
Z(@)= 2, Aue 
To derive the relationships between the measured time series and the unknown 


frequency-directional wave spectrum the cross spectrum must be considered. In general, 
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the cross spectrum between two time series X(t) and Y(t) with Fourier transforms X(@) 


and Y(q@) is defined as 


(X (@)Y (w)) 


Cy (@) = Wy 


(7.11) 


where " indicates the complex conjugate, [Soong 1993]. Substitution of Equation 7.10 
into Equation 7.11 yields 
C yy (@) = >) cos @ sin 6E(@, 6) (7.12) 
6 


where it is assumed that the wave components propagating in different directions are 


Statistically independent. The cross spectrum Cyy can be expressed in continuous form as 


2h 
Cy (@) = | cos@ sin 6E(@, 0)d0 . (7.13) 


0 
Cross-spectra of the other time series pairs can be obtained in a similar manner. The full 
set of relations for buoy cross-spectra can be found in [Dean 1984]. 
It is convenient to define a normalized directional distribution of energy at a 
frequency Was 
_ E(@,@) 


S(@; ; 
(0; @) E(o) (7.14) 
with unit integral 
Lr 
se | E(@,0)0 ae 
SCI a 7.15 
| E(@) E(@) —_ 


With this definition, Equation 7.13 and the other referenced spectral relations can be 


combined and expressed in terms of $(@;q@). Dropping the frequency dependence these 


relations can be expresses as 


ime) _ on = 
(Cy +Cy Cy) = cos Sad ad, 


Im(C,, ) ae 
(Gate, cay Bown 
) ia : Estee zi (7.16) 
es = | cos265(@)d0 =a, 
C 2c, ; 
2x 

sh = | sin 205(0)d@ = b, 

XX YY. 0 


These four relations between the cross-spectra of the buoy measurements and the 
directional distribution of wave energy, derived by [Long 1980], form the basis for most 


of the buoy analysis techniques currently used. 


D. EXTENSION TO SUBSURFACE SENSORS 


As discussed in the previous section, the motion of a wave buoy 1s directly related 
to the fluid velocity, therefore, the cross-spectra of a tri-directional current meter yields 
the same low resolution directional wave information obtained from buoy measurements. 
Substituting the normalized spectra of the vertical (Z) velocity component w, and the 
horizontal (X, Y) velocity components u and v into Equation 7.16, the lowest four Fourier 


moments of the directional distribution of wave energy can be obtained and are given by 


sO eee 
\O- occ me 
a,(Q@) oe ans (7.19) 
b,(Q) a (7.20) 


where C(q@) is the spectral matrix of the velocity components u, v, w. Since the wave 
direction, @ , is referenced to the navigation frame (N-E-D), vehicle borne sensor 


measurements must be transformed prior to use. It is interesting to note that the estimates 
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of these directional moments are insensitive to errors, so long as the errors are the same 
on all measurement axes of the sensors, which is typical with oceanographic sensors 
installed on AUVs. 

The objective of the data analysis is to infer the directional distribution S(@), from 
the four measured moments a), b;, a2, and bo. The most widely used techniques are 


described below. 
a) The Cosine Power Distribution 


[Longuet-Higgins 1963] introduced a simple cosine-power distribution, 
S(0) = Acos”* (FF 2p 


with @nean the mean propagation direction, s a parameter that controls the width of the 
distribution and A, a normalization coefficient. The parameters Qnean and s are 
determined by fitting Equation 7.21 to the relations given in Equations 7.17-7.20. 

The main drawback to this simple method is that Equation 7.21, with only two 
free parameters can describe only unimodal distributions, and thus fails in situations with 
a bimodal sea state (e.g., multi-directional seas during a wind veering event or swell 


arriving from two different sources). 
b) The Maximum Entropy Method 


[Lygre and Krogstad 1986] introduced the maximum entropy or MEM 
estimate of S(@). Unlike Equation 7.21, this approach can describe both unimodal and 
bimodal distributions and exactly fits the relations given in Equations 7.17-7.20. This 


directional spectrum is given by 


i= * = * : 
5(8) = ee (7.22) 
As (I-g,e — 9, 


with 
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Gent, 


C, =a, +ib, 


, wa Suen (7.23) 
Mears 2 

1-|c,| 
Q, =C,—C,9, 


Stull, the directional distribution is poorly constrained by only four moments and the 


estimates require careful interpretation, [Krogstad 1991]. 
c) Mean Direction and Directional Spread 


An alternative approach that avoids the pitfalls of S(@) estimation, is to 


describe the directionality of waves by a few simple parameters. For narrow S(@), a 
mean propagation direction @, and a root-mean-square measure of the directional 
spreading energy O, can be defined in terms of the first-order and second-order Fourier 


moments a,, b,, a, and b, [Kuik 1988], given by, 


ie tls 2. C21) 
a, 
ao; =2{1-[a, cos(O, ) +b, sin(@, )]] 22) 
ate [2 (7.23) 
2 a, 
OC; al - la, cos(20, )+b, sin(2@, )|]. (7.24) 


Again, this method fails to identify bimodal distributions but it is useful to 
determine a base direction so that a control command could be determined. More on this 


approach will be discussed later in this chapter. 


E. | CORRECTION FOR A MOVING PLATFORM (THE DOPPLER 
EQUATION) 


The equations for the wave directionality estimations presented in the previous 
section is valid for a non-moving sensor. However, when information is obtained from 


an AUV, corrections must be made to account for the vehicle motion. As discussed in 
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Chapter II, the wave frequency which the vehicle encounters while moving through the 
wave field has been shifted. This shift can be determined by using the well known 
Doppler equation, Equations 2.56 or 4.25. 

Using the techniques outlined in the previous section will give the wave 
directional distribution as a function of vehicle encounter frequency. If these estimations 
are used in conjunction with the Doppler equation in a recursive manner, the estimation 
of S(@) can be found as a function of true frequency. 

| The method by which this is performed is outlined below; 


e Determine the three components of the fluid velocity in vehicle body fixed 
coordinates. } 


e Transform the fluid velocities into the global navigation frame using Equation 
Zeon 


e Compute the auto- and cross-spectra of the fluid velocity components. 
e Determine the Fourier moments using Equations 7.17-7.20. 


e Convert the Fourier moments into Krogstad notation and use the MEM to 
determine the directional distribution S(@). 


e Using the vehicles mean velocity, and the mean direction obtained from 
Equation 7.21 or 7.23, apply the Doppler equation to determine the frequency 
shift. 

e Return to 3 and complete the process until frequencies converge. 

The corrections to the estimation of $(@) using the Doppler equation are quite 
small for slow moving vehicles. Considering, for example, the NPS Phoenix AUV which 
has a maximum forward velocity of 1.5 m/s, the error associated with not using the 
Doppler equation are approximately + | sec in identification of wave periods between 4 
and 40 seconds. Similarly, the error in direction estimation is approximately 5-7 degrees. 
When an AUV goes into a station keeping mode, where vehicle velocities are 


significantly reduced, the modifications required due to the Doppler shift are negligible. 
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F. DETERMINATION OF CONTROL COMMANDS 


The resulting directions and variances, from Equations 7.21-7.24, are for each 
frequency component of the wave field. To use this information in determining a heading 
command, single direction must be found. Bulk Fourier moments, weighted by the 


energy density of the wave field, 


a’ a, (O) 
b? 1 Fig b, (@) 
emer) 7.24 
Cae J J a,(Q) a 
b; b,(@) 
with E” the swell variance given by, 
oe 
E’ = | E(aydo, (7.2#) 
fi 


may be substituted into Equations 7.21-7.24 to yield a bulk fluid direction and variance. 
It is this bulk fluid direction that is used in generation the heading command dependent 


on the mission requirements. 


G. SUMMARY 


This chapter has presented the techniques currently employed for the 
determination of wave directional estimations from standard wave following buoys. It 
has extended this analysis for use in determining directional estimates from data gathered 
by an Autonomous Underwater Vehicle. Using this data gathered, an approach was 
presented which will allow the deployed vehicle to obtain information about the 
directionality of its working environment thereby allowing the vehicle to have 


information available to make decisions regarding mission execution. 
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VHl. EXPERIMENTAL RESULTS 


A. INTRODUCTION 


With the theoretical presentation of this dissertation complete, the primary focus 
of this chapter will be on the experimental validation of the Disturbance Compensation 
Controller (DCC). In addition, results from the ONR sponsored AUVFEST 798 will 
prove the concept of wave direction estimation from a moving platform presented in 
Chapter VII. 

This chapter will begin with a discussion of the real-time implementation of the 
DCC in the NPS Phoenix AUV. It will follow with a presentation of the experimental 
results, conducted in Monterey Harbor, which displays the achievable performance of the 
DCC. Comparison between theoretical and experimental results will be made. Lastly, 
directional energy and spectral plots obtained by Phoenix, while operating in the Gulf of 


Mexico, south of Gulfport, MS, will be shown. 


B. SOFTWARE IMPLEMENTATION 


The implementation of this control process is unique since it is split between the 
two CPUs installed in Phoenix. The NPS AUV uses a Pentium based PC-104 running 
QNX and a GESPAC Card Cage running OS9 for mission control and execution. The 
DCC requires information from both processors, connected by Ethernet sockets, to 
compute and pass the commanded propeller RPMs to the execution level. 

The control architecture presently running in Phoenix is based on shared memory 
processes. The PC-104 computer runs a “main” process that controls the synchronization 
of the data sharing, while the GESPAC clock controls the real-time contro] features. The 
two-processors use the shared memory as the common data buffer, controlled by 
semaphores to ensure the information transfer is consistent with the clock speed. A 
graphical representation of this description is shown in Figure 8.1. As seen in the 
graphic, for the DCC implementation, all needed processes are run in the PC-104 with the 
main purpose of the GESPAC is to provide the commanded voltages to the propulsion 


motors. 
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Figure 8.1 Software Implementation of DCC 


A block diagram of the DCC implementation in the Phoenix AUV was given in 
Chapter VI, Figure 6.1. It is redisplayed as Figure 8.1 for easy reference. This diagram 
represents the melding of the software and the hardware in the vehicle. The ground 


velocity is from the RDI, the relative velocity from the ADV and y, r from the 


directional gyro. 
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Figure 8.2 Block Diagram of DCC Implementation 
Sy 


c. EXPERIMENTAL VALIDATION OF THE DCC 


The DCC was tested in Monterey Harbor between the months of March and May 
1999. During this time, the Phoenix was held under surge control for over 90 minutes, 
during various runs, without a drive off. Table 8.1 provides a sample of the runs 


conducted -—- the validation of the controller. 





Date rh aad or) ae Pore 
4/2/99 4 min 0.3624 2.96 high gain, 
Pel eat eee ae 
0.6324 high gain, 
short run, 


vehicle physical 


disturbed 
0.4312 0.265 high gain, 

he 
0.5090 0.285 high gain, 


hal oe is min 0.5508 0.108 high gain, single 
shaft 


ie. a min 0.3620 0.192 medium-high gain, 
ADV noise 
problem 
ft ll 
meal min 0. - 0. aed low gain, ADV 
_ ae 
10 min 0.3587 0.202 medium-high gain, 
ADV noise 
Ceeer = 
10 min 0.4276 medium-low gain, 
ADV noise 
problem 


Table 8.1 Sample Summary of DCC Validation Runs 













medium-low gain, 
ADV noise 






problem 
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Defining a measure of performance, the disturbance rejection ratio (DRR), as the 
ratio of standard deviation of the vehicle ground velocity to the standard deviation of the 
fluid velocity, the ability of an AUV to reject disturbances for different conditions and 
control designs can be compared. Referring to the DRR definition, for perfect 
disturbance cancellation the DRR will be equal to zero, while for designs where the 
control input excites the vehicle, as was shown in Chapter V, the DRR will be greater 
than one. For each operating point, the standard deviation of the propeller response is 
normalized by the maximum propeller revolutions, max. 

Table 8.1 indicates that excellent disturbance rejection was achieved, even for the 
short runs where only limited statistical information was recorded. The tests showed that 
even when the vehicle was disturbed by a source other than the fluid velocity, it was able 
to return to the commanded position in a stable fashion. 

A series a plots, Figures 8.3-8.8, show the results of one of the validation runs. 
This run was conducted on April 22, 1999 in Monterey Harbor. The Phoenix was 
commanded to a navigational position of 0 meters in the longitudinal direction. As the 


results indicate, the vehicle behaved as expected. The standard deviation of the 
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Figure 8.3 Comparison of Measured and Estimated Position, April 22, 1999, Run*3 
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Figure 8.5 Measured Ground Velocity, April 22, 1999, Run*3 
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Figure 8.6 Comparison of Measured and Estimated Relative Velocity, 


Run*3 
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Figure 8.7 Fluid Velocity Estimate, April 22, 1999, Run*3 
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Figure 8.8 Estimated Thrust, April 22, 1999, Run’3 


positional error was 9.6 cm with ground velocity standard deviation of 1.5 cm/s. 

This run was the most interesting of the validation runs conducted. At the 
beginning of this run, it was noticed that the starboard shaft was not turning. Even with 
this propulsion system casualty, the vehicle was able to hold position and the controller 
did not go unstable. With only one shaft turning the effective input gain for the vehicle 
was reduced by 50%. Operations of this nature indicate a very robust design. It can be 
seen in Figure 8.4, that there is a small increase in propeller revolutions around the 50 
second point of the run. Data analysis indicated that this was approximately when the 
starboard shaft failed. Investigation into the cause of the shaft failure determined that a 
universal joint in that shafting had worked loose. 

As a graphical representation of the performance expected for the DCC a various 
simulations ware conducted, using the asynchronous simulator discussed in the Chapter 
VI, with the estimated fluid velocity obtained during this run (April 22, 1999, run” 3) as 
the disturbance. The gains on the DCC were varied to produce a position response verses 
propeller response curve. The actual experimental results, presented in Table 8.1, were 


superimposed on the theoretical curve obtained from simulation. These results are shown 
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in Figure 8.9. As seen, the experimental and theoretical results are very close indicating a 
physically realistic simulator. 

The comparisons displayed in Figure 8.9, contrast data analyzed from 
experimental results together with computer predicted results with the same disturbance, 
and yield insight into the achievable performance of the DCC. It indicates that there is a 
limit to the amount of disturbance rejection that is physically realizable. This limit is 
controlled by ability of the propulsion system to produce the needed input to maintain 
position without excessive oscillations. The excessive oscillations have a detrimental 
effect of the life of the propulsion system. 

As a note, the short runs, displayed in Figure 8.9, were conducted with a 
controller gain parameter of 7 = 100, a high gain. If the length of these runs were 


extended, these points would move closer to the curve as with the other runs displayed. 
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Figure 8.9 Comparison of DCC Performance, Simulation and Experimental 
Up to this point, the only results presented are for the Phoenix maintaining 
position to the origin, the point which the run was initiated. Questions arise as to how 


effective the controller is in dealing with transients. This question may be answered by 


referring to Figure 8.10. This figure depicts the transient response of the Phoenix for the 
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various DCC gains presented in Figure 8.9. As the figure indicates, the controller deals 
with transients extremely well. 

The responses displayed in Figure 8.10 are for the regulator solution. What is 
meant by this, recalling that the SMC formulation requires kinematically consistent 
position, velocity and acceleration, is that no command inputs, other that position were 
used. In doing this, it is expected that the vehicle will over shoot and oscillate around the 


commanded position consistent with some settling time. 
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Figure 8.10 Comparison of Transient Response for Various Control Gains 


With these transient responses come some issues with regard to operational 
implementation. If the goal is to position the vehicle close to, but with out touching, an 
object, some means of predicting the transient must be available. A resulting question 
that needs to be answered is; Does the developed simulator, which, based on the 
comparison in Figure 8.9, accurately predict the transient response? By comparing the 
results of the experimental runs and the simulated results, for the same disturbance input 
and DCC design, see Figure 8.11, the question can be answered, "yes.” 

As seen in this plot, the simulated results accurately reflects the measured 
transient response of the Phoenix. The steady state response, however, does not match 


exactly. The reason for this is two-fold. First, the Phoenix, for recovery reasons, is 


is? 


maintained approximately two-pounds buoyant. This weight and buoyancy mismatch, as 
discussed in Chapter II, cause additional excitation forces resulting from the wave 
induced fluid accelerations. Since the fluid acceleration cannot be measured, this 
additional excitation force is difficult to replicate in simulation yielding errors between 
the real and simulated response. Second, the experimental results are measured from a 
6DOF rigid body, where as the simulator results come from a 1DOF surge model. The 


coupling effects from the surge-pitch dynamics will effect the comparison. 
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Figure 8.11 Transient Response Prediction of the DCC 


D. WAVE DIRECTIONAL ESTIMATION USING THE NPS PHOENIX AUV 


During November 1998, the NPS Center for AUV Research, under the direction 
of Professor Anthony Healey, participated in the ONR sponsored AUVFEST '98. This 
AUV technology demonstrator was held in the Gulf of Mexico, south of Gulfport, MS. A 
complete description of the exercise can be found in Appendix D. 

The Phoenix AUV was used during this exercise as a mobile sensor to gather 
oceanographic data. Using the concepts presented in Chapter VII, the vehicle conducted 
27 short-term sampling missions. The products that were obtained included directional 


energy plots, directional spectrum plots and mean current estimations. 
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The key to the ability of Phoenix to obtain data capable of producing this 
information is the combined ADV/RDI/MotionPak/Directional Gyro sensor suite. With 
these sensors, accurate three dimensional fluid velocities, expressed in global quantities, 
were capable of being obtained in post-processing. Since the RDI/ADV sensors are 
collocated, little vehicle induced motion remains after processing the data. 

The data obtained validated the concept of obtaining tactical oceanographic data 
from an underwater vehicle. During the collection of the data, remnants of Hurricane 
Mitch were still present in the Gulf, providing an excellent source of ground swell. In 
addition, there was a significant wind generated wave component in a different direction 
than the swell component, resulting in a multi-modal spectrum. 

The results presented in Figures 8.12-8.14 provide a sample of the oceanographic data 
obtained during this offshore exercise. As seen in Figure 8.12, the bi-modal properties of 
the seaway are captured, as well as an estimate of the significant wave height. The ability 
to estimate the dual directions is due to the use of the MEM algorithm presented in 
Chapter VII. Figure 8.13 presents the associated spectrum plots for this energy estimate. 
The ability of an AUV to estimate currents is shown in Figure 8.14. Using a triangular, 
time based run, the current can be determined using set and drift calculations from the 
error in final position as well as the heading error on each leg. This information can be 
feed directly into the vehicle’s navigation process to account for the offset due to current 
thereby increasing navigation accuracy. Short term averages from each of the three legs 


are in general agreement with the overall average determined from the navigational drift. 


E. SUMMARY 


This chapter has presented the experimental results associated with this 
dissertation. The results validate the theory and show that it is possible for a small AUV 
to reject wave induced disturbances making it capable of performing positioning tasks in 
shallow water. In addition, the robustness of the design to sensor noise and propulsion 


faults has been demonstrated. 
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Phoenix Survey Data 
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Figure 8.12 Sample Direction Energy Plot From Phoenix Data, Nov. 4, 1998 (Run” 9), 
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Figure 8.13 Sample Direction Spectrum Plots From Phoenix Data, Nov. 4, 1998 
(Run” 9), Gulf of Mexico 
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As a supplement, due to the sensor suite available, it was shown that tactical 
oceanographic data is obtainable from a moving AUV. The vehicle in this manner 
becomes an intelligent, mobile off-board sensor, thus presenting the argument for AUV 


deployment with operational fleet units. 
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Figure 8.14 Short-term Current Estimation from Phoenix, November 8, 1999, (Run’ 2), 
Gulf of Mexico | 
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IX. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 


This research has shown through modeling, simulation and experimental 
validation that intervention tasks performed by intelligent underwater robots are 
improved by their ability to gather, learn and use information about their working 
environment. The development, implementation and verification of a real-time 
embedded Disturbance Compensation Controller (DCC) for small AUVs, has provided a 
new technology showing that it is possible to use underwater vehicles for station-keeping 
tasks in shallow water. 

The work conducted in support of the dissertation objectives has provided three 
specific contributions to the field of tetherless underwater robotics. First, a new 
generalized approach to the modeling of small underwater vehicles subject to shallow 
water wave and current effects was developed. Using appropriate modeling 
formulations, the fluid disturbance effects are incorporated directly into the equations of 
motion leading to model based control laws for disturbance compensation. In addition, 
this formulation provides the ability to construct a generalized distributed simulation 
capability for the evaluation of underwater vehicle systems in shallow water, which is 
particularly useful to U.S. Navy tactical decision making. 

Secondly, it is proven using asynchronous simulations and in-ocean experimental 
validations, that substantial compensation of wave induced disturbances may be achieved 
from direct on-line measurements of the water column velocities. This technique 
eliminates the need to develop and incorporate sophisticated predictive disturbance 
models in the control system design. 

Thirdly, it is shown how small underwater vehicles, using direct fluid 
measurements, can obtain short-term wave magnitude and direction, as well as current 
estimates, thereby providing useful information in the area of tactical oceanography. It is 
also shown how a general seaway direction may be obtained from this information for 


use in mission planning and control. 
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In addition to the general contributions listed above, several specific conclusions 


may be drawn from this research. In particular, 


The input requirements associated with the DCC are vehicle position, relative 
velocity, propulsion force, and water column velocity. The attractive nature 
of the DCC is that these quantities may be measured or estimated from a hull 
mounted sensor suite. With these the full benefits of DCC can be realized. 


The DCC proved to be very robust to sensor noise and propulsion system 
faults. Stable vehicle performance was maintained in the presence of the loss 
of a propulsion shaft. 


The design of the propulsion system must allow for a rapid, oscillatory 
response to the call for propeller speed commands. In addition, considerations 
for fault tolerant operations are recommended. 


For the application of dynamic station-keeping, propulsion system lags and 
associated thrust reduction terms must be identified and taken into account. 


Identified parameters of the nonlinear model provided stable and easily 
tunable controllers, as verified by in-water experiments. Excellent agreement 
was achieved between experimental and simulated results. 


Seaway models developed using AR representations require high order to 
effectively achieve spectral matching. However, lower order models (four or 
six) can be used adequately for estimation and control. 


Extended Kalman filtering methods for seaway estimation appear to be 
satisfactory. 


RECOMMENDATIONS FOR FUTURE RESEARCH 


As a result of the work performed in this dissertation several research areas have 


arisen requiring further investigation. These include: 


The validation of the 6DOF model for other control modes with field data is 
recommended. This is particularly interesting in flight control where motion 
minimization Is critical to improve side-scan imagery. 


Application of the disturbance compensation techniques presented herein to 
other control modes used in shallow water positioning needs inquiry. By 
achieving compensation in all three planes, small AUVs may be employed in 
a number of positioning tasks, including mine neutralization. 
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A sensitivity study as to the required degree and accuracy of the disturbance 
model, with regard to disturbance compensation, is warranted. It is believed 
that improved disturbance dynamics models may increase DCC robustness. 


While this work has used EKF methods for the identification of seaway 
dynamics, other techniques, such as neural networks and MUSIC, may have 
advantages that could be explored. 


The DCC formulation does allow for prediction of future water column 
velocities. This added information may possibly provide additional benefits 
not explored here, although, this study has indicated that the DCC is not 
highly sensitive to future information. 
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APPENDIX A. EQUATIONS OF MOTION AND PARAMETERS FOR THE NPS 
PHOENIX 


The equations of motion, parameter description and parameter values used to 


simulate the dynamic behavior of the NPS Phoenix AUV is given in this Appendix. 


PHYSICAL PARAMETERS 


[Faramser [Bessie [va 
Weight 1934.9 N 

a cae aime 
mass 197.2 kg 
Buoyancy 1934.9N 

eee Length C2 ih 









Mass Moment of Inertia 3.66 N-m-sec 
about x-ax1s Cy ft-Ib-sec”) 

Mass Moment of Inertia 56.94 N-m-sec 
~ about y-axis (42. ft-Ib-sec’) 


Mass Moment of Inertia 61.01 N-m-sec* 
about z-axis (42.0 ft-Ib-sec’) 


oes Cross Product of Inertia 0.0 N-m-sec 
ae about xy-axes (0.0 ft-Ib-sec’) 
Ly, Cross Product of Inertia 0.0 N-m-sec” 
nea about yz-axes (0.0 ft-Ib-sec’) 






Le Cross Product of Inertia 0.0 N-m-sec 
about xz-axes (0.0 ft-Ib-sec*) 
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x Coordinate of CG From 
Body Fixed Origin 

y Coordinate of CG From 
Body Fixed Origin 

z Coordinate of CG From 
Body Fixed Origin 

x Coordinate of CB From 
Body Fixed Origin 

y Coordinate of CB From 


Body Fixed Origin 
z Coordinate of CB From 
Location of Bow Vertical 
Location of Stern Vertical 
Location of Bow Lateral 
Thruster from CG 


Location of Stern Lateral 
Location of Left Propeller 


Location of Left Propeller 
from CG 
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0.003 m 
(0.01 ft) 
0.003 m 
(0.01 ft) 
0.0128 m 
(0.042 ft) 
0.003 m 
(0.01 ft) 
0.00 m 
(0.0 ft) 
0.00 m 
(0.0 ft) 
0.432 m 
(1.420 ft) 
-0.432 m 
(-1.420 ft) 
0.585 m 
(1.920 ft) 
-0.585 m 
(-1.920 ft) 
-0.10 m 
(-0.33 ft) 
0.10 m 
(OS3rt0) 


CONTROL INPUTS 


a ei. 
(+ 5 Ibs) 









(+ 10 lbs) 
+44.45N 
(+ 10 Ibs) 







Right Propeller Force 
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(+ 5 lbs) 
te ON 


(+ 5 Ibs) 








Force 








Bow Vertical Thruster 









Force 





Bow Vertical Thruster - to 25 IN 













Force 








Stern Plane Deflection 
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NON-DIMENSIONAL HYDRODYNAMIC COEFFICIENTS 
Surge Hydrodynamic Coefficients: 


76 













rears | serena | ees | 


Sway Hydrodynamic Coefficients: 
Y a =0.01241 


Y;=0.0 Y, =-0.03430 
Y/ =-0.00178 Y’=0.0 Y’ =0.0 Y,,=0.01241 


a 


Heave Hydrodynamic Coefficients: 


Z;, =-0.00253 Z’, =-0.09340 Z’ =-0.78440 
dp 


Y,,, =0.0 


















Roll Hydrodynamic Coefficients: 
K, 














K,=0 (0 









2 


Pitch Hydrodynamic Coefficients: 


M ; =-0.00625 M‘,=-0.00253 M’ =0.05122 


Yaw Hydrodynamic Coefficients: 


==U)- = : : 
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EQUATIONS OF MOTION 


The equations of motion for the NPS Phoenix AUV are presented in this section 





of the Appendix. These equations are based on the vector-matrix equation given in 
equation 2.87. The expanded equations use the assumption that Phoenix possesses both 


horizontal and vertical plane symmetry. 


Mass Matrix: 
m—X, OF 0 0 MG —My>¢ 
0 m-Y, 0 ~(mzg +Y,) 0 mx, -Y, 
re 0 0 im (be, MY¢ (mx, +Z,} 0 
| 0 =(mzg + K;) myc I,,-K, -I,, —(I,,+K;) 
m2, 0 —(mx, + M,) -I,, I, -M, re 
hie mx, -N, 0 -(,,+N,) 9 -I, ie 
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Centrifugal/Coriolis Matrix: 





C(x)= 
[ 0 — mr mg m(¥oq + Zr) —migqg-Z.w —migr+Y.v 
mr 0 — mp — mygp+Z.w m( 2¢r + X¢ Pp) -myor—X ju 
mq mp 0 —micp-Y,v —miggt X ju m(x¢ Pp + eq) 
= m(yoqQ+ Zr) myo p-Z.w migpr+Y,v 0 —1.q-loptac aN Lr t+Typ-lyqt+M oq 
micqt+Z,w —-m(Z6r + Xp) mzoqg—-X,u Ld tle — ieee 0 Hr -lyq+l,p-K;p 
mxor-Y,v my or + X ju —M(X¢P + ¥oq) atta ttg—-M 7g fol g—i,p+ Kp 0 
Control] Allocation Matrix: 
B(x)= 
UrX 55, +UU]X 5,5, 55, UGX asp, +UUU|X syysyy Sop UX rs, +UlulXs,5,5,, UgXq5,, tuXe6 5, 1 1 #0 0 0 0 
ululYs, 0 ululY.,, 0 oo Oo 1 OR 
0 ujulZ», 0 uulZ.,,, 0.0 1.0” 7a 
0 0 0) 0 0 0 0 0 oO oO 
0 ululM », 0 ulu|M 0 0 =x, 0 =x mmo 
ulu|V ., 0 ululN ,, 0 -y, -y, O x, O x, 


Control] Input Vector: 


u=[0,,,0p,0 0 Leer a Tye roll 


sr? vi? 


Sp? 


Euler Angle Rates and Global Positions: 
X= U, +ucosycosé@ + vicosysin Osin @— sin ycos@]+ w[cos ysin @cos@ + sinysing] 
Y= V,+u sin ycos@ + v[sin ysin @sin ¢ +cosycos@¢] + usin psin@cos¢ — cos ysing] 


Z =W, —usin@+vcos@sin ¢ + wcos@cos¢ 


Q= p+qsin@tanG@ + rcos gtanO 
0=qcosg-—rsing 


.  (gsin@+rcos@) 
a 
cos@ 
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Rotation Matrix: 


This rotation is based on equation 2.86. 


gion. 2 
10 = T, 


cosycosé sin wc@ 
T, =| cosysin @sin g—sin ycos@ sin ysin Osin @+cosycos@ 
cosysin@cos¢@+snysng  smysin@cos@—cosysingd 


l 0 —sin 6 
T,= 0 cosé sin @ 
0 —sng  cos@cos¢ 


Surge Motion Equation: 


m[u — vr +wgq —x,(q° + r* )+ VAG = 1) be a) eg oe 


=D +X q° +X r° gfe Nepean ar Ve Ger 0 te ce ete Nea 
+ U(X 95, Or + X 45, Os, ) + U(X 15, Op, + Xi5,9,) 


2 2 2 2 2 2 
+X,v° +X,,,w tulul(Xs5 (6, +6, eo + X55 Or, ) 


-~(W -—B)sin0@+F,,+F,,-X.,.u 


res 





ul) 


W -B \c- ) 
(7 hee, tT yV, +2 3W, +10, +2V, +7,W,| 
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~—sin@ 
cos @sin @ 
cos@cos @ 


Sway Motion Equation: 


m[V+ur—wp+Xg( pqt+r)—yol p +r )+2%6(qr- p)]+Z,wp +X ,ur 
=V pri ty pat), q +Y iv) up ty unt Va kW ae 


+Y,uv+Y,,vwtululY,,6, +¥s,6,, ) 
£ Coll x Nv + ar vt xr)| dx 


+(W -B )cosOsin@+F,, +f ,, 


(= ew, aaa +T,W, +TU , +TrV, + TW, | 


Heave Motion Equation: 


m[w—uq+vp+x¢( pr- g)+yolart+p)—-2%6( Pp’ + q° )]-Y,vp + X ,ur 
ee: 2 2 : 
=Z6G + Zp mace pt + Lee ate. Wed een ale Ce 


Zep + Z,,uw + ulul( Z sp Op Cr, ) 
= [Cab x ow — 2x9 ow - xq) am 


TY ot COSOCOS a) mete 


(=? le, +Ty.V,+TyW, +150 5 + TV + TW, | 
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Roll Motion Equation: 


I,pt+d, - 1,)ar+1,(pr - Q-1,lq° -r°)-1,(pqt+ r)-Z,vw+Y,vw 
—N,qr+M ,qr+mly.(w—uq+ vp) — Z¢(v + ur — wp)] 


cae ee Dd akg 
+K,v+K,up+ Kur+K,vqg+K,,wp+K,,wr+ K,uv+ K,,vw 


+(y,W — y,B)cosécos@ — (z,W — z,B)cos@sin @ 
Wz, -B é : , 
“| ig ~PXp Jew, +TyV,+Ty3W, +T,U ,+TyV, +T,,W, | 
§ 


Wy. —By, \- : . 
[Pom iu ptTyV,+TW, +TyU ,+ToV, +T WV, | 
g 


Pitch Motion Equation: 


. . ° 2 2 
Lg +, 1, pr! (Ges py tA a) te 
—m[x¢(w—uq+vp)—z,4(u—vr+wq)]+Z,uw—-X ,uwt+N, pr—K, pr 
aS = 2 2 . 
=M,q+M,,p +M,pr+M_r°+M,w+Muq+M,vp+M,,vr 


+ M,v’+M,uw+ulul(M 6, +M 5,6, ) 


C,.b¢ x w— xq w - xq )x]dx 


Sve = hp) COSO COs =a WV — ces) SING aaah). as ane 
Wz. —Bz : ; : . ; . 
-[e= Fee Iu Pie eta) eee Ore eee tan 


Wx,.—-Bx, \r- : 
(Ro eu, TLV + Lag W yp + 15) U y +LaV, + TW, | 


Lee 


Yaw Motion Equation: 


Lr+(1,-1,)pa-1,(p° —@° )—1,(pr+q)+1,(ar—p) 

+ m[x,(v + ur — wp) y,(u—vr+wq)] —Y,uv + X,uv-M,pq+K,pq 
=N,pt+N,r+N,,pqt+N,,gr+Nv+N,up+N,ur+N,vqat+N, wp +N,,wr 
+ N,uv + N,,vw +ulul N;,6,, + N3,5,, ) 


. [C,,h¢ xv + xr Nv + x7 )x] dx 


tail 


+(x,W —x,B)singcos@ + (y,W — y,B)sin@ 
Mp ho a” ae oe Vee eg 


Wyv.-B : 
[o-Ps Mu p+TyV, +TW, +70, +T2V, +T,W, | 


(Ae — Bx, 


Jw ptTV, +TyW,+TU ,+TV, +T,W, | 
g 
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APPENDIX B. DOPPLER SENSORS 


This Appendix provides an overview of the two Doppler sensors used for control 
implementation in this dissertation; namely the SonTek® ADV Ocean and the RDI® 


Navigator DVL. 


THE SONTEK ADVOCEAN 

The SonTek ADVOcean (Acoustic Doppler Velocimeter Ocean Probe) is a 
versatile, high-precision instrument used to measure 3D water velocity in the most 
challenging applications, Figure B.1. The ADVOcean is designed for a wide range of 


environments including the surf zone, open ocean, rivers, lakes, and estuaries. 
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Figure B.1 ADVOcean Probe 


The ADVOcean uses acoustic Doppler technology to measure 3D flow in a small 
sampling volume located a fixed distance (18-cm) from the probe, Figure B.2. The 
velocity range is programmable from +5 to +500 cm/s. Data can be acquired at sampling 
rates up to 25 Hz. 

With no zero offset, the ADVOcean can be used to measure flow velocities from 
less than 1 mm/s to over 5 m/s. The remote sampling volume, stability, and rapid 
response of the ADVOcean make it perfect for al 1 types of flow measurement: mean 


currents, waves, and turbulent flow parameters. 
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The ADVOcean consists of two elements: a probe and processor. The probe 
includes the acoustic sensors, receiver, and optional sensors in a submersible housing. It 


is connected to the processor using a custom shielded cable. 






Remote Sampling ‘ 
Volume for 
3D Velocity 


Figure B.2 ADV Sampling Volume 


The ADVOcean processor operates from external DC power and outputs data 
using serial communication or a set of analog voltages. The processor can be operated 
from any PC-compatible computer or can be integrated with a variety of data acquisition 
systems. 

STANDARD ADVOCEAN 

The standard ADVOcean probe, Figure B.3, is designed for long-term 
deployments in hostile environments. The entire probe is constructed from 316 stainless 
steel, and protected from corrosion by a sacrificial zinc anode. With no moving parts, the 
ADVOcean has excellent resistance to biological fouling. For added protection, the entire 
probe, including the transducers, can be coated with anti-fouling paint. The probe is 


connected to the processor using an underwater mateable connector. 
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Figure B.3 Standard ADVOcean Probe 


For deep-water deployments, the ADVOcean can be rated for depths up to 2000 
m (the standard depth rating is 400 m). Deep-water ADVOcean systems are commonly 
used on autonomous underwater vehicles (AUVs) and remotely operated vehicles 
(ROVSs) for detailed microstructure measurements. 

In any configuration, the ADVOcean probe is immune to zero drift and has no 
inherent minimum detectable velocity. The probe calibration can only change with 
physical damage to the system. No regular maintenance or re-calibration 1s needed. 
ADVOCEAN WITH OPTIONAL SENSORS 

The ADV Ocean probe can include a number of optional sensors to greatly expand 
its Measurement capabilities. These include a compass and 2-axis tilt sensor allowing the 
ADVOcean to report velocity data in Earth (East-North-Up) coordinates; a pressure 
sensor for wave height estimation (PUVW) and surface-level measurements; and a 
temperature sensor for automatic sound speed compensation. 

ADVOcean probes with optional sensors use an expanded housing, Figure B.4, 
constructed from acetyl (Delrin), and have the same reliability and performance as the 


standard ADVOcean probe. 
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Figure B.4 ADVOcean Probe with Optional Sensors Housing 


ADVOCEAN PROCESSOR 

The ADVOcean processor can be housed in two different ways depending upon 
whether the processor will need to be submerged. The ADVOcean processor operates 
from DC power and is typically connected to a portable computer running SonTek's 
powerful data acquisition software. It can also be integrated with a variety of data 
acquisition systems using serial communication or the analog output voltages. 
STANDARD FEATURES 


ADVOcean systems include the following standard features listed in Table B.1. 


ADVOcean Probe ADVOcean Processor 


Factory Calibration (can only change Dual serial communication (RS-232 
standard, RS-422 for cable lengths to 
1500 m) 















with physical damage to the probe) 





e Programmable velocity range from 


+5 to +500 cm/s 
Submersible to 400 m 








e Four analog output voltages (3 








velocity, 1 signal strength) for 














e 10-m cable to processor (50-m integration with analog data 












max.) acquisition systems 






e Hardware synchronization with 





external sensors using sync input and 






output 





Table B.1 Standard ADVOcean Features 
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OPTIONS 

Several options are available for use with ADVOcean systems. The most 
common are the compass, pressure and temperature sensors. The compass and 2-axis tilt 
sensors allow the ADVOcean to report velocity data in Earth (East-North-Up) 
coordinates. The sensor has a built-in calibration feature to compensate for magnetic 
distortion. The user can easily re-orient the compass for up, down, or side-looking 
operation. Any ADVOcean with compass/tilt or pressure sensor includes a temperature 
sensor to compensate for changes in sound speed. Sound speed is used to convert 
Doppler-shift to water velocity. 
ADVOCEAN PERFORMANCE SPECIFICATIONS 

The performance specifications of the ADVOcean are listed in Table B.2. More 


information may be found at the SonTek web site http://www.sontek.com. 


Performance Specifications 


General Operation Compass/Tilt Sensor 
Acoustic frequency: 5 MHz e Resolution (heading, pitch, roll): 0.1° 
Sampling rate: Programmable from 0.1 to | ¢ Accuracy (heading): +2° 
25 Hz : Accuracy (pitch, roll): +1° 
Sampling volume size: 2.0 cm3 Temperature Sensor 
Distance to sampling volume: 18 cm @) Resolution. jMice 


Minimum water depth: 20 cm © Accuracy: +0.1°C 


Input power : 12-24 VDC Pressure Sensor 


Velocity Data e Available full-scale ranges: 10, 20, 50, 
Range: Programmable to +5, 20, 50, 200 100, 200, and 500 m 


or 500 cm/s ¢ — Resolution: 0.00025 x (full scale) 


Resolution: 0.01 cm/s e Accuracy: +0.5 percent full scale 
Accuracy: +1 percent of measured | Environmental 


velocity, +0.25 cm/s e Operating temperature: -5° to 40°C 


Random noise: Approximately 1% of | , Storage temperature: -10° to 50°C 


velocity range at 25 Hz; Dimensions 


@ See Figure B.5 





Table B.2, ADVOcean Performance Specifications 


je 


42.3 cm / 
16.7 in 







Diameter 
1.3 cm/ 
0.5in 


rere eee tee ee ee ee ee ee ee ee ee 


a 
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Diameter.” 
20.1 cm/ 
7.9 in? 


Figure B.5 ADVOcean Dimensions 


THE RDI NAVIGATOR DVL 

The WORKHORSE NAVIGATOR DVL is designed to provide rapid, precise 
velocity updates. Its small size, high accuracy, and low power consumption make it well 
suited to applications such as station keeping and sea bed surveys from underwater 
vehicles. The NAVIGATOR can be integrated with existing navigational systems 
(USBL, LBL and/or inertial). Its high-resolution velocity data provides a better focused 
picture of your vehicle's location and altitude. The NAVIGATOR is less than half the 
length and weight of standard broadband DVLs. Average power consumption of the 
WHN-1200 is only 8 watts. The NAVIGATOR is about 60 % of the price of standard 
broadband DVLs. The NAVIGATOR uses the patented Broadband technology. The 
WHN-1200 kHz has velocity accuracy of +/-0.2% +/- 0.2 cm/s. The NAVIGATOR 
measures Velocity, altitude, heading, pitch/roll, and temperature. For a description of the 


operations of this sensor refer to [Gordon 1996] 
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WHN-1200 SPECIFICATIONS 
The specification for the 1.2 MHz Navigator are presented in Table B.3. Further 
specifications can be found at the RDI web site http://www.rdinstruments.com. 


Transducer and Pressure Case 


Actual Frequency 1229 kHz 
Beamwidth eZ 
Beam Angle (from vertical) 30° 
Configuration 4-beam-convex 
Housing & Transducer 6061 aluminum 
Material 
External Connector 7-pin low-profile 
Weight (in air) 6.4 kg 
Weight (in water) 4.2kg 
Altitude 
Minimum 03 m 
Maximum 30 m 
Bottom Velocity 
Short Term Error 0.3 cm/s 
(V = 1.0 m/s) 
Long term Error +0.2% +0.2 cm/s 
Ping Rate 1-10 Hz 


Table B.3 Navigator DVL Specifications 
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APPENDIX C. THE NPS PHOENIX AUV 


The physical layout of the NPS Phoenix AUV is shown in Figure C.1. Detailed 
description of the vehicle can be found in [Marco 1996] and [Brutzman 1997]. In 
addition an online description can be found on the NPS center for AUV Research web 


Site at http://www.cs.nps.navy.mil/research/auv/auvstats. 





ST725 SONAR 
RDI DOPPLER ST1000 SONAR 
SONAR 
SonTek ADV 
DEPTH CELL 
TRANSDUCER 
FIN SERVO (8) 
BOW LEAK 
DETECTOR BOW LATERAL 
ADV PROCESSOR THRUSTER 
BOW VERTICAL MOTION PAK 
THRUSTER 
DC POWER SUPPLIE 
COMPUTER POWER AND bans : 
SUPPLY (2) 
MOTOR SERVO 12 VOLT BATTERY (2) 
CONTROLLERS (6) FOR COMPUTERS 
FreeWave RADIO 
COMM. LINK GESPAC COMPUTER 
WITH SIGNAL BOARDS 
QNX PENTIUM 
COMPUTER 
12 VOLT BATTERY (2) 
STERN VERTICAL FOR GYROS/MOTORS 
THRUSTER 
FREE GYRO EL SHe 
STERN LATERAL DETECTOR 
THRUSTER 


REAR SCREW MOTOR (2) 


CONTROL FINS (8) REAR SCREW (2) 


Drawn by D. Marco 1999 
Figure C.1 Physical Layout of Phoenix AUV 


Prior to September 1996, Phoenix was used extensively as a test tank research 
vehicle. To transition the vehicle and the center to an ocean going operational unit 
required some significant upgrades in vehicle capabilities and center acquisitions. Table 


C.1 summarizes the upgrades and acquisitions that were performed to allow Phoenix to 
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become deployable in the ocean. It also highlights the design features of the new vehicle, 


Figures C.la and C.1b, presently being outfitted at the center. 





Figure C.la Wire Frame Diagram of New NPS AUV 





Figure C.1b Solid Model of New NPS AUV 
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Pre-July 1996 Present New Vehicle 


Propulsion Motors two 24v brushed DC two 1/4 hp, 24v brushless two 48v DC, 1/4 hp, 
DC (see figure C.2) brushless 


Propellers 3.5 in model submarine 7 in ducted propellers (in 7 1n ducted propellers (in 


















house designed) (see figure house designed) 


propellers (see figure C.3) 
. C.4) 







Control Actuators seven control surfaces, (one seven control surfaces, (one 









eight control surfaces, (two 













fwd rudders, two aft rudders, lower fwd rudders, two aft lower fwd rudders, two aft 


and a pair of bow and stern rudders, and a pair of bow rudders, and a pair of bow 


















planes), four 3 .5 in thrusters and stern planes), four 3 .5 and stern planes), four 6 in 


(two vertical fore and aft and in thrusters (see figure thrusters (two vertical fore 





two horizontal fore and aft) C.4a(two vertical fore and and aft and two horizontal 










aft and two horizontal fore fore and aft) 
and aft) (see figure C.5) 
GESPAC Computer System | PC-104 with a Pentum chip 
(166 MHz) running QNX 


(see figure C.6) 


(see figure C.5) 











Vehicle Control GESPAC Computer System 


running OS-9 real-tume 















Computer running OS-9 real-ume 







operating system operating system 





PC-104 with a Pentum chip | PC-104 witha Pentum chip 


(166 MHz) running QNX 


Mission Control Sun Solaris 








(90 MHz) running QNX (see 






Computer 
figure C.6) 


Electrical Power two independent 24v lead two independent 24v lead single 48v Absorbed Glass 
Ballast System manual lead ballast manual lead ballast with Variable ballast system 
i i ee eee 


Communication Ethernet cable while in test Ethernet cable while in test Ethernet cable while in test 
System tank tank, 900 MHz spread tank, 900 MHz spread 


(see figure C.6) 





spectrum modem dunng spectrum modem and u/w 


missions (see figure C.7) acoustics during missions 





Navigation System DR using water speed and EKF fusing Doppler and EKF fusing Doppler, GPS, 


vehicle heading vehicle motion LBL and vehicle motion 
Attitude Sensor three axis Mechanical Rate 6DOF Solid State inertial 6DOF Solid State inertial 
Gyro's and a vertical gyro sensing system (see figure sensing system 
C.8) 


Table C.1 Vehicle Upgrades and Acquisitions 
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Heading Reference Directional Gyro Directional Gyro with a Directional Gyro (possible a 


PrecisionNav electronic Honeywell ring laser gyro) 
compass backup (see figure with a PrecisionNav 
C.9) electronic compass backup 


Speed Reference turbo probe RDI DVL and SonTek ADV {| RDI DVL and SonTek ADV 
(see figures C.10 and C.11) 


Sensor Suite ST-725 Scanning Sonar, ST-725 Scanning Sonar, ST-725 Scanning Sonar, 
ST-1000 Imaging Sonar (see ST-1000 Imaging Sonar, ST-1000 Imaging Sonar, 
figure C.12) RDI DVL, SonTek ADV RDI DVL, SonTek ADV, 


















Altimeter , video camera 


Support equipment Transport cart (see figure Wells Cargo® Travel trailer | Wells Cargo® Travel trailer 


C.13) outfitted with a mobile lab, outfitted with a mobile lab, 
12' inflatable boat and 12’ inflatable boat and 
trailer, portable generator trailer, portable generator 
and computer workstation, and computer workstation, 


vehicle shipping containers vehicle shipping containers 
LXT acoustic tracking 
system (see figures C.14, 


C15, C6 andealy) 


Table C.1 Vehicle Upgrades and Acquisitions (Cont.) 
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Figure C.2 Brushless DC Motors 





Figure C.3 Old 3-in. Props 


eh 





Figure C.4 Present 7-in Ducted Props 
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Figure C.4a Present 3.5-in Thrusters 
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Figure C.5 Phoenix Control Actuators 
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Figure C.6 Mission Control Computer 
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Figure C.7 FreeWave Modem 





Figure C.8 Systron Donner MotionPak 
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Figure C.9 PrecisionNav TCM2 Compass 





Figure C.10 RDI Navigator DVL 
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Figure C.12 ST-725 and ST-1000 Sonars 
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Figure C.13 Transport Cart 
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Figure C.14_ Mobile Lab Internals 
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Figure C.16 Support Craft 


io 


S 


ere RAS 


EWA ORSAY 





Aa 


* 
Pe 


os 


“3 er teapmeceninirs tiled 





Figure C.17 LXT Acoustic Tracking System 
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APPENDIX D. AUVFEST '98 


OVERVIEW 

The Naval Postgraduate School (NPS) participated in the NAVO sponsored 
AUVFEEST for the first time, [Bunce 1998]. The lessons learned from this exercise were 
extremely valuable for future operational planning and vehicle development. First, the 
NPS Phoenix vehicle, a research platform, had never been deployed offshore from a 
research ship before. This task provided challenges that the four-man team failed to 
recognize beforehand, but was able to overcome. Like all other participants, the vehicle 
and support equipment had to be transported from their base of operations to Gulfport, 
MS. This movement of equipment was new to the Center for AUV Research 
(www.cs.nps.navy.mil/research/auv), but enlightened the Center on the logistics involved 
with offshore operations. 

During the work-ups for this exercise, two significant hardware problems 
prem First, the Doppler unit originally integrated into Phoenix failed the week just 
prior to departure for Gulfport. This sensor was beyond repair, and a new Doppler 
needed to be purchased. A RDI Navigator DVL, a $25K unit, was purchased, through 
the Naval Supply System, and delivered to the Center in five days. The purchase of this 
unit in this time period was remarkable considering the government regulations that must 
be followed for a major purchase of this type. The vehicle, the support equipment, the 
new Doppler and all other sensors were shipped to Gulfport, installed, integrated into the 
vehicle control system and tested without faults in three days, a major accomplishment 
for a group that had never performed missions away from their base of operation. 

Secondly, during the weekend prior to shipboard load out, after the vehicle had 
been reassembled and all systems had been successfully tested and verified as 
performing, the SonTek Acoustic Doppler Velocimeter (ADV) was physically damaged. 
The damage to the ADV was not recognized until the after the vehicle had been loaded 
aboard the R/V GURE, and vessel was underway. The Center was able to contact the 
vendor and have a new unit shipped overnight to Gulfport where immediate configuration 
and installation was accomplished without error, affording Phoenix the ability to remain 
operational. 
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The objectives of the Center during AUVFEST 798 were two-fold. First, NPS 
planned to use the sensors suite installed on Phoenix to demonstrate the ability of a 
moving platform to characterize the seaway. The theory and results from this section 
have been presented in Chapter VII and VIII, respectively. The second goal was to use 
the Phoenix’s forward-looking, sector-scanning sonar (Tritech ST725), to image water 
column targets and demonstrate the Phoenix’s onboard target identification algorithm. 
Navigation of Phoenix was accomplished using the RDI DVL together with a directional 
Syro heading reference in a deadreckoning filter. Details with regard to the above goals 
' and vehicle missions are presented in the following sections. 

LOGISTICS 

The Phoenix AUV and all its associated support equipment was air freighted from 
Monterey, CA to Gulfport, MS via FedEx. This was a challenge for the Center since this 
was the first operational deployment. Custom shipping containers were purchased from 
Hardigg Industries, see Figure D.1, to hold the hull and nose fairing. These containers 
performed extremely well and protected the vehicle. The NPS shipping department 
provided packaging of the centers support equipment, see Figure D.2. Again due to the 
professional nature of the NPS employees, not a single item was damaged or tured up 


missing. 





Figure D.1 Phoenix Shipping Containers 
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Figure D.2 Packaging of Support Equipment 


Upon arrival in Gulfport, the Center set-up shop in the NAVO BOATDET office, 
see Figure D.3. The vehicle was reassembled, including the newly purchased Doppler, 


bench tested and water tested, without any system degradations in two days. 





Figure D.3 Gulfport Temporary Work Space 


SHIPBOARD OPERATIONS 

| The Phoenix, with its support equipment, was loaded aboard the Texas A&M 
University research ship, R/V GYRE, see Figure D.4. In addition to NPS, there were 
participants from Florida Atlantic University (FAU), Woods Hole Oceanographic 
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Institute (WHOI) and Massachusetts Institute of Technology (MIT) taking part in the 


demonstrations. 





Figure D.4 R/V Gyre at Anchor 


During the seven day exercise, the Phoenix conducted 27 separate open ocean 
runs, see Table D.1 for a sample. This was a significant accomplishment seeing that the 
Phoenix was designed for tank testing some 10 years earlier. 

Operations from the ship were challenging. This was the first time Center 
personnel launched Phoenix from a crane, see Figure D.5. By the end of the second day 


of operations, the Phoenix crew was acting as if this was “old hat.” 





Figure D.5 Phoenix Being Launched from the Gyre 


The exercise was conducted in the Gulf of Mexico, south of Gulfport, MS. Figure 
D.6 displays the operating area. In this operating area, there were three different mission 
boxes. Each box had separate features so that the vehicles from the various participants 


could demonstrate their capabilities. 
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Figure D.6 AUVFEST Operations Area 


DIRECTIONAL WAVE INFORMATION 

Directional wave information was obtained during the AUVFEST missions. See 
Chapters VII and VIII for more information. 
AUTOMATIC SONAR TARGET RECOGNITION EXPERIMENTS 

One of the objectives of the AUVFEST experiments involved testing a recently 
developed sonar target recognition software module. The software is designed to process 
the sonar returns in real time and provide a reduced data set from the large amounts of 
data gathered. From the data, the centers of concentrations of high intensity sonar returns 
are identified and their locations stored during mission time, which requires no post 
processing. The location information can be then be used for post mission analysis or 
path re-planning to return to these areas for further study during the same mission. High 
intensity concentrations suggest underwater objects or structures, while areas of low 
intensity do not. Since the majority of the open ocean is clear, only a small amount of the 
data gathered is meaningful, and this approach greatly reduces the data storage 
requirements of the onboard computer system. The following will give a brief description 
of the sonar used, the identification algorithm, and finally the results from one of the 


missions. 
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Day/Data | Run comments on Run Problems/ Fixes Sea State comments General info 
a eee name | Number 


salinity 29.5 ppt at surface to 
34.0 ppt at bottom. Water 
temp ~74 degrees. 




















added buoyancy in nose / 
tow float/reset fins/ 


light (1 foot seas, multi- 
directional, wind ~10 mph) 


1. ‘Set up/balancing/venicle 
was heavy by nose, ~2lbs 
light. Time base run 4 mins 

(2 min out, 2 in back) at 3’ 

depth. 


































sea state 1-2 (2-3 feet, 4 foot 
swells). Seas from multi- 
directions. Swells from 115 
true. Wind driven from 025 
true. Wind 10-15 mph 








Taking ships heading/go to 
Gyre doasmall WP run, a 
large WP run, then a run to 
area 2 buoy 


TB run. At 090 true to go to) Fix sign error in DG offset 
GB. Initialize dg to ships _| calc. Due to sign error, we 
head. Went to command _ {ended up with twice the 
but 180 off expected ship's heading, which was 
‘course. Mag compass 330 | 100 degrees, error. 
- ; degrees. 
1104.03 2 repeat of first run. no 
changes to script. Same 
’ results as first run. 
1104.04 | 3 . Same runas first. 
‘Changed DG offset by 
2*pi. Results were no 
| different than first run 
Same run as first. 
Attempted to give a 270 
heading command. Results 
were no different than first 
run 
‘Reinitialized vehicle Behaved as commanded. 
headed at Gyre. New zero |Run was successful. 
for DG. (~085 true) Time 
based run for 2 minutes 
with a commanded 
peading of 000. 


















{Fix sign error in DG offset 
calc 













Fix sign error in DG offset 
calc 

















Fix sign error in DG offset 
calc 

















1104_06 















Checked battery (24.5v) and 
computer (22.8v) voltages 

when mission was complete. 
Leak detectors at 1.09. 


Behaved as commanded. Comments from the chase 
Run was successful. boat was that we were doing 
~2.5 moh. 





five minutes. Heading Run was successful. 
‘away from Gyre. 






= Time based run at -090 for |Behaved as eominanticd 



















Time based run at 045 for 
seven minutes. Heading 
‘towards Gyre. 
Time based run at 045 for 
| two minutes. Heading 
towards the Gyre. 


1104_09 








Behaved as commanded. 
Run was successful. 












seconds heading away 
from the Gyre. Way points 
were (0,40,30), (-40,80,3), 
(-70,110,0) relative to 
initialized heading 
long waypoint run. 
Attempted to start at (0,0) 
and returned to (0,0). 
Vehicle timed out trying to 
| get to (100,-85) due to 
turning the long way. 































Waypoints were (0,0), (50,0), 
(100,0), (150, -10), (190,-40), 
(200, -60), (180,-80), (140,- 
90), (100,-85), (50,-65), (20,- 
30), (0,0). With 6-meter 
diameter watch circle. 


Run was moderately 
successful. 




























Waypoints were (0,0), (50,0), 
(100,0), (150, -10), (190,-40), 
(200, -60), (180,-80), (140,- 
90), (100,-85), (50,-65), (20,- 
30), (0,0). With 10-meter 
diameter watch circle 


Run was moderately 
| Attempted to start at (0,0) j|successful. Fix is to correct 
and returned to (0,0). the heading command with 
Vehicle timed out trying to |an if statement to ensure 
get to (100,-85) due to that the vehicle always 
turning the long way. turns the shortest direction. 
Attempted to increase 
watch circle. 


long waypoint run. 













Table D.1 Sample Phoenix Missions 
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SONAR HEAD GENERAL DESCRIPTION 

The NPS Phoenix is equipped with two mechanically scannable high frequency 
sonar heads built by Tritech Corp. One is a ST725 scanning sonar and the other is a 
ST1000 profiler sonar. Each head can be scanned continuously through 360 degrees of 
rotation or swept through any defined angular sector around the central axis of the unit. 


During normal operation, the head will ping, wait for return echoes to process, and then 


rotate a specified angular width and repeat. Step widths of 0.9 i 1.8, and 3.6 are 
computer selectable. 
All missions performed at AUVFEST °98 used the ST725 which operates at a 


frequency of 725 kilohertz and emits an acoustic beam 2.5 * wide in the horizontal plane 


by 24° wide in the vertical plane. This device is described as a scanning sonar due to the 
nature of the range and intensity information returned for each ping. A scanning sonar 
operates by placing the intensities of the echoes from each ping into discrete segments of 
range called range bins. For this sonar, the number of range bins is nominally 128, but for 
operating ranges of 10 meters or less, the number of range bins is reduced to 64. The 
maximum operating range of the ST725 is 100 meters with a minimum operating range 
of 6 meters, and provides a resolution of (Maximum Range)/128 or (Maximum 
Range)/64, depending on the range setting used. The range setting used in the Gulf was 
20 meters, which gave a resolution of approximately 15 cm. 

In order to more clearly analyze the returns, the data can be thresholded to 
analyze only returns above a certain intensity level so that significant objects/structures 
can be seen, while other less significant entities (e. g. suspended particles in the water 
column, weak multi-path echoes, noise, etc.), are excluded. Combining thresholding with 


an appropriate power setting for the transducer, high quality results can be achieved. 
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SONAR CLUSTER IDENTIFICATION PROGRAM 
The identification algorithm is designed to recognize areas of contiguous high 
intensity sonar returns. Below 1s a section of a test case showing sonar scanlines that 


have been thresholded to record intensities above 10. 


0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000011110000000000000000000000 

000000011 13121100000000000000000000 
000000011 13 15 13 12 1100000000000000000 
00000012110 141312 1211000000000000000 
000000012 13 141413 12 121100000000011000 
0000000000000000000000000000000000 
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00000000 12 1314 130000000000000000000 

0000000001100000000000000000000000 

0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000014000000000000000000000000000 
0000001300000000000000000000000000 

0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
0000000000000000000000000000000000 
00000000132110000000000000000000000 
00000000 12 13 14 1314130000000000000000 
000000013 13115 15 1112 1212 1200000000000 
000000001213 141314000000000000000000 
00000000001111000000000000000000000 
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000000000 13130000000000000000000000 
00000000141400000000000000000000000 
00000001111000000000000000000000000 
00000000121200000000000000000000000 
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The algorithm records the centroid (X, Y pairs) of each cluster of high intensity 
returns, while ignoring noise or small concentrations such as the 14, 13 group shown 
above. Several parameters are selectable to tune the algorithm for target identification 
such as maximum cluster width, breadth, number of non-contiguous contacts, etc. The 
following presents the results of the identification algorithm from a run at AUVFEST. 
SONAR RESULTS FROM AUVFEST 

Since there were no submerged targets in the area where the Phoenix operated, the 
chase boat served as a Suitable target. Figure D.7 shows the targets identified (chase 
boat) with the centroid of each marked with a cross-hair. The trajectory of the Phoenix is 
shown with a solid line while the sonar returns with an intensity above 10 are shown with 
asterisks. Figure D.8 shows the lower right target cluster with the centroid clearly 
identified by the algorithm. Since the AUVFEST results were very positive, further 
experiments will be conducted in Monterey Bay in the early summer of 1999. 
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Figure D.7 Clusters Identified. 
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